*COMDECK BASE
      COMMON/BASE/NMAX,JMAX,KMAX,LMAX,JM,KM,LM,DT,GAMMA,GAMI,SMU,FSMACH
     1  ,DX1,DY1,DZ1,ND,ND2,FV(5),FD(5),HD,ALP,GD,OMEGA,HDX,HDY,HDZ
     2,RM,CNBR,PI,ITR,INVISC,LAMIN,NP,INT1,INT2,INT3
     3,KPLANE,ITW,TW,WMDOT,L1BC,LMAXBC,LU,LL,JWL1,KU,WTFAC,JMAXBC
     4,J1BC,JB,KMAXBC,KAL,K1BC,
     5DTFAC,RESID,NK,AREF,LW,KW,JL1WL,JL1WU,JK1WL,JK1WU,JLW,JKW,INXBC
     6,DQ0QMX,JU,RMACH,DTMAX,KVIS,LVIS,KLVIS
     7,HTOT(500),PTOT(500),VOU(500),WOU(500)
      COMMON/GEO/NB1,NB2,RFRONT,RMAX,XR,XMAX,DRAD,DXC
      COMMON/READ/NRST,IWRIT,NGRI
      COMMON/VIS/RE,PR
*COMDECK BTRID
      COMMON/BTRID/A(60,5,5),B(60,5,5),C(60,5,5),D(60,5,5),F(60,5)
*COMDECK COUNT
      COMMON/COUNT/NC,NC1,TAU,NROUT
*COMDECK RHSBCC
      COMMON/RHSBCC/IW,KUORLU,LORKMX,KORLMX,LMORKM,LLORKL,LUORKU,LKMXBC
     1,KLMXBC,LK1BC,KL1BC,LORKW,KORLW,JLK1WL,JLK1WU,JKL1WL,JKL1WU,JLKW
     1,JKLW,DZORDY,GKPR,PRTR,I,J,KLORLL
*COMDECK VARS
      COMMON/VARS/Q(500,6,60)
      COMMON/VAR0/S(500,5,60)
      COMMON/VAR1/X(500,60),Y(500,60),Z(500,60)
      COMMON/VAR3/XX(60,4),YY(60,4),ZZ(60,4)
*COMDECK VISCO
      COMMON/VISCO/RMUEXP,ISUTH,TSUTH
*DECK DRIVER
      PROGRAM NOZL3D(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,
     1 TAPE2=ICFILE,TAPE4=NOZOUT,TAPE3=NOZIN)
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
      DIMENSION GF(5)
C
C
C         CORE STORAGE VERSION OF 3D NOZZLE FLOW PROGRAM
C        3-D UNSTEADY PARABOLIZED NAVIER-STOKES EQUATIONS SOLVED WITH
C        IMPLICIT ALGORITHM
C         THIS IS THE MAIN CONTROL ROUTINE
C
C
C
C  INITIALIZATION
C
      NC1 = 0
      TAU = 0.
      RESIDS=0.0
      NRES=1
      NK=0
      CALL INITIA
      IT0 = NC1
      ITMAX = NC1 + NMAX
      KL0 = (LMAX-1)*ND
C
C  COMPUTE MAXIMUM EIGENVALUE AND COURANT NUMBER
C
      CALL EIGEN
      REIN=RE*FSMACH
C
C         WRITE OUT INPUT DATA
C
      WRITE(3)KMAX,JMAX,LMAX,ITMAX,LMAXBC,L1BC,FSMACH,GAMMA,REIN,SMU,DT,
     1                                                       ALP,CNBR,PR
      WRITE(3) IT0,TAU,((X(K,J),Y(K,J),Z(K,J),K = 1,KMAX),J = 1,JMAX)
      WRITE(3) IT0,TAU,((X(KL0+K,J),Y(KL0+K,J),Z(KL0+K,J),K = 1,KMAX),
     1   J = 1,JMAX)
C
C         COMPUTE GLOBL PERFORMANCE PARAMETERS USING INITIAL CONDITIONS
C
      CALL GLOBL(GF)
      WRITE (6,10)IT0,TAU,DT
   10 FORMAT(5X,5H NC =,I4,7H TIME =,F11.4,5H DT =,E12.4)
      WRITE(3) IT0,TAU,GF
      WRITE(4)KMAX,JMAX,LMAX,ITMAX,LMAXBC,L1BC,FSMACH,GAMMA,REIN,SMU,DT,
     1                                                       ALP,CNBR,PR
C
C         OUTPUT INITIAL FLOWFIELD
C
      IF(IWRIT.EQ.1) CALL OUTALL(NC1)
      CONTINUE
C     MAIN LOOP MAIN LOOP MAIN LOOP MAIN LOOP MAIN LOOP MAIN LOOP
      DO 60 N=1,NMAX
      NC = N+NC1
C
C  IMPLICIT INTEGRATION
C
      CALL STEP
C
C  CALCULATE AND WRITE CURRENT GLOBL PERFORMANCE PARAMETERS
C
      CALL GLOBL(GF)
      WRITE(3) NC,TAU,GF
C
C  OUTPUT VARIABLES EVERY NP TIME STEPS
C
      IF (N-(N/NP)*NP) 35,35,50
   35 CONTINUE
      WRITE (6,40)NC
   40 FORMAT(1H0,5H*****,3HN =,I5,7H  *****)
      IK=(KMAX+1)/2
      CALL OUT2(1,1)
      CALL OUT2(1,LMAX)
   50 CONTINUE
      NK=NK+1
C
C         OUTPUT FLOWFIELD EVERY NROUT TIME STEPS
C
      IF (MOD(N,NROUT).NE.0) GO TO 55
      IF (N.EQ.NMAX) GO TO 60
      CALL OUTALL(NC)
   55 CONTINUE
      IF (N.EQ.NMAX) GO TO 60
C         BUMP TIME STEP,DT,AND SMOOTHING COEFFICIENT SMU BY THE FACTOR
C         DTFAC WHENEVER DQ0QMX IS LESS THAN OMEGA AND DT IS LESS THAN
C         DTMAX.  SMU IS BOUNDED ABOVE BY 0.8 FOR STABILITY.
C
      IF(ABS(DQ0QMX).GT.OMEGA) GO TO 60
      IF (DTFAC.EQ.1.) GO TO 60
          IF(DT .GE. DTMAX) GO TO 60
      DTF=DTMAX/DT
      IF(DTFAC.GT.DTF) DTFAC=DTF
      IF(ABS(DTFAC-1.0) .LT.1.E-4) DTFAC=1.0
      DT=DT*DTFAC
      SMU=SMU*DTFAC
      IF(SMU.GT.0.8) SMU=.8
      HD = .5*DT
      HDX = HD/DX1
      HDY = HD/DY1
      HDZ = HD/DZ1
      CNBR=0.0
      CALL EIGEN
   60 CONTINUE
      NC1 = NC
      CALL OUTALL(-NC1)
      STOP
      END
*DECK AMATRX
      SUBROUTINE AMATRX(A,J,KL,R1,R2,R3,R4)
*CALL BASE
*CALL VARS
      DIMENSION A(5,5)
      GAM2 = 2.-GAMMA
      RR = 1./Q(KL,1,J)
      U = Q(KL,2,J)*RR
      V = Q(KL,3,J)*RR
      W = Q(KL,4,J)*RR
      UU = U*R1+V*R2+W*R3
      UT = U**2+V**2+W**2
      C1 = GAMI*UT*.5
      C2 = Q(KL,5,J)*RR*GAMMA
      A(1,1) = R4
      A(1,2) = R1
      A(1,3) = R2
      A(1,4) = R3
      A(1,5) = 0.
      A(2,1) = R1*C1-U*UU
      A(2,2) = R4+UU+R1*GAM2*U
      A(2,3) = -R1*GAMI*V+R2*U
      A(2,4) = -R1*GAMI*W+R3*U
      A(2,5) = R1*GAMI
      A(3,1) = R2*C1-V*UU
      A(3,2) = R1*V-R2*GAMI*U
      A(3,3) = R4+UU+R2*GAM2*V
      A(3,4) = -R2*GAMI*W+R3*V
      A(3,5) = R2*GAMI
      A(4,1) = R3*C1-W*UU
      A(4,2) = R1*W-R3*GAMI*U
      A(4,3) = R2*W-R3*GAMI*V
      A(4,4) = R4+UU+R3*GAM2*W
      A(4,5) = R3*GAMI
      A(5,1) = (-C2+2.*C1)*UU
      A(5,2) = R1*(C2-C1)-GAMI*U*UU
      A(5,3) = R2*(C2-C1)-GAMI*V*UU
      A(5,4) = R3*(C2-C1)-GAMI*W*UU
      A(5,5) = R4+GAMMA*UU
      RETURN
      END
*DECK BC
      SUBROUTINE BC
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL RHSBCC
      COMMON/AXISYM/LAXIS
      DIMENSION RUVWE(2,5),PP(2),T(2),V(2),W(2)
C     DIMENSION A1(50),B1(50),C1(50),D1(50),F1(50)
      P(I,J)=GD*(Q(I,5,J)-.5*(Q(I,2,J)**2+Q(I,3,J)**2+Q(I,4,J)**2)
     1   /Q(I,1,J))
      EAVG(IL,IB,J)=.5*(Q(IL,5,J)+Q(IB,5,J))
      EORAVG(IL,IB,J)=.5*(Q(IL,5,J)/Q(IL,1,J)+Q(IB,5,J)/Q(IB,1,J))
C
C ELIMINATE ROUNDOFF ERROR IN TRANSVERSE VELOCITY V FOR PLANE FLOW
C     IN THE XZ PLANE(KPLANE=1) AND ENFORCE PLANAR SYMMETRY.
C
      IF (KMAX.EQ.1) GO TO 95
      IF (KPLANE.EQ.0) GO TO 95
      DO 90 J=1,JMAX
      DO 90 L=1,LMAX
      KLM=(L-1)*ND
      DO 90 K=2,KMAX
      KL=KLM+K
      KL2=KLM+1
      DO 85 N=1,5
   85 Q(KL,N,J)=Q(KL2,N,J)
      Q(KL,3,J)=0.0
   90 CONTINUE
   95 CONTINUE
C     LAGGED B.C. FOR INTERNAL CORNERS AND SINGULARITIES
C     INTERNAL CORNER AT KMAX,LMAX
      IF (KPLANE.EQ.1) GO TO 105
      IF (KMAXBC.NE.1.OR.LMAXBC.NE.1) GO TO 105
      KL=(LMAX-1)*ND+KMAX
      KLL=KL-1
      KLB=KL-ND
      DO 100 J=JB,JMAX
      Q(KL,5,J)=EAVG(KLL,KLB,J)
  100 Q(KL,1,J)=Q(KL,5,J)/EORAVG(KLL,KLB,J)
  105 CONTINUE
      IF (NC.EQ.0) GO TO 140
C
C     ELIMINATE ANY OUTFLOW AT COMPUTED INFLOW PLANE J=1
C
C     IF (J1BC.NE.1) GO TO 115
C     DO 110 L=LL,LU
C     DO 110 K=KAL,KU
C     KL=(L-1)*ND+K
C     IF (Q(KL,2,1).GT.0.) GO TO 110
C     Q(KL,2,1)=0.
C     Q(KL,3,1)=0.
C     Q(KL,4,1)=0.
C     Q(KL,5,1)=PTOT(KL)/GD
C     Q(KL,1,1)=PTOT(KL)/HTOT(KL)
C 110 CONTINUE
C 115 CONTINUE
C      B.C AT SINGULAR AXIS L = 1 FOR AXISYMMETRIC FLOW (LAXIS=1)
      IF (LAXIS.EQ.0.OR.KPLANE.EQ.1) GO TO 140
C     SET V=W=0 AND EXTRAPOLATE U,RHO,T QUADRATICALLY WITH ZERO NORMAL
C     GRADIENT AT AXIS
      DO 135 J=JB,JMAX
      RHO=0.
      U=0.
      TEM=0.
      DO 120 K=1,KMAX
      KL=K
      KL1=K+ND
      KL2=KL1+ND
      RHO1=Q(KL1,1,J)
      RHO2=Q(KL2,1,J)
      U1=Q(KL1,2,J)/RHO1
      U2=Q(KL2,2,J)/RHO2
      V1=Q(KL1,3,J)/RHO1
      V2=Q(KL2,3,J)/RHO2
      W1=Q(KL1,4,J)/RHO1
      W2=Q(KL2,4,J)/RHO2
      T1=GD*(Q(KL1,5,J)/RHO1-.5*(U1**2+V1**2+W1**2))
      T2=GD*(Q(KL2,5,J)/RHO2-.5*(U2**2+V2**2+W2**2))
      RHO=RHO+4.*RHO1-RHO2
      U=U+4.*U1-U2
  120 TEM=TEM+4.*T1-T2
      U=U/(3.*KMAX)
      RHO=RHO/(3.*KMAX)
      TEM=TEM/(3.*KMAX)
      IF (J.NE.1) GO TO 125
      IF (J1BC.EQ.0) GO TO 125
      TEM=HTOT(1)*(RHO*HTOT(1)/PTOT(1))**GAMI
      U=SQRT(2.*(HTOT(1)-TEM)/GAMI)
  125 CONTINUE
      DO 130 K=1,KMAX
      Q(K,1,J)=RHO
      Q(K,2,J)=RHO*U
      Q(K,3,J)=0.
      Q(K,4,J)=0.
  130 Q(K,5,J)=RHO*(TEM/GD+.5*U**2)
  135 CONTINUE
C     ELIMINATE ROUNDOFF ERROR FOR AXISYMMETRIC FLOW
C     PI2=2.*ATAN(1.)
C     PIP=PI2/(KMAX-1.)
C     DO 600 J=JB,JMAX
C     DO 600 L=1,LMAX
C     KLO=(L-1)*ND
C     KL1=KLO+1
C     DO 599 K=2,KMAX
C     KL=KLO+K
C     THET=(K-1)*PIP
C     Q(KL,1,J)=Q(KL1,1,J)
C     Q(KL,2,J)=Q(KL1,2,J)
C     Q(KL,5,J)=Q(KL1,5,J)
C     Q(KL,3,J)=Q(KL1,4,J)*SIN(THET)
C 599 Q(KL,4,J)=Q(KL1,4,J)*COS(THET)
C     Q(KLO+KMAX,4,J)=0.
C 600 CONTINUE
  140 CONTINUE
C     LAGGED B.C. FOR SINGULAR CORNERS AND EDGES
C     AXIAL INTERSECTION OF VERTICAL AND HORIZONTAL WALLS IN PLANES
C     L=1 AND K=1
      IF (JL1WU.LE.0.OR.JK1WU.LE.0) GO TO 150
      J1=MAX0(JL1WL,JK1WL)
      J2=MIN0(JL1WU,JK1WU)
      KL=1
      KLL=KL+1
      KLB=KL+ND
      DO 145 J=J1,J2
      Q(KL,5,J)=EAVG(KLL,KLB,J)
  145 Q(KL,1,J)=Q(KL,5,J)/EORAVG(KLL,KLB,J)
C     CORNERS AND EDGES OF INTERMEDIATE WALLS
  150 IF (KW.LE.0.AND.LW.LE.0) GO TO 185
      IF (KW.GT.0.AND.LW.GT.0) GO TO 170
C     TRAILING EDGE OF LONE INTERNAL WALL THAT SPANS ENTIRE WIDTH OF
C     COMPUTATIONAL DOMAIN
      DO 165 IW=1,2
      CALL SW
      IF (LORKW.LE.0) GO TO 165
      J=JLKW
      I1=0
      LORK=LORKW
      ISHIFT=-ND
      IF(IW.EQ.2) ISHIFT=-1
  155 DO 160 K=KLORLL,KUORLU
      KL=(LORK-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+LORK
      KLB=KL+ISHIFT
      Q(KL,1,J)=Q(KLB,1,J)
  160 Q(KL,5,J)=P(KLB,J)/GD
      IF (I1.EQ.1) GO TO 165
      I1=1
      LORK=LORKW+1
      ISHIFT=-ISHIFT
      GO TO 155
  165 CONTINUE
      GO TO 255
C     CORNERS AND EDGES OF BOTH INTERMEDIATE WALLS
  170 CONTINUE
C     AXIAL INTERSECTIONS OF THE TWO WALLS
      J2=MIN0(JKW,JLW)
C     INTERNAL CORNER
      KL=(LW-1)*ND+KW
      I1=0
  175 KLL=KL-1
      KLB=KL-ND
      IF(NC.EQ.NRST) GO TO 177
      Q(KL,5,1)=PTOT(KL)/GD
      Q(KL,1,1)=PTOT(KL)/HTOT(KL)
  177 CONTINUE
      DO 180 J=2,J2
      Q(KL,5,J)=EAVG(KLL,KLB,J)
  180 Q(KL,1,J)=Q(KL,5,J)/EORAVG(KLL,KLB,J)
      IF (I1.EQ.1) GO TO 185
C      EXTERNAL CORNER
      I1=1
      KL=KL+ND+1
      GO TO 175
  185 DO 250 IW=1,2
      CALL SW
C     INTERSECTION OF PLUG AT LORK=1 WITH SIDEPLATE AT KORL=KORLW OR KORLMX
      IF (JLK1WU.LE.0) GO TO 200
      LORK=1
      KORL=KORLW
      IF(KORLW.EQ.0)KORL=KORLMX
      IF(KORL.EQ.1) GO TO 200
      KL=(LORK-1)*ND+KORL
      KLB=KL-1
      KLL=KL+ND
      IF (IW.EQ.1) GO TO 190
      KL=(KORL-1)*ND+LORK
      KLB=KL-ND
      KLL=KL+1
  190 DO 195 J=JLK1WL,JLK1WU
      Q(KL,5,J)=EAVG(KLL,KLB,J)
  195 Q(KL,1,J)=Q(KL,5,J)/EORAVG(KLL,KLB,J)
C     TRAILING EDGES OF WALL AT J=JLKW
  200 IF(JLKW.EQ.0) GO TO 250
      J=JLKW
      K2=KORLW-1
      LORK=LORKW
      I1=0
      ISHIFT=-ND
      IF(IW.EQ.2) ISHIFT=-1
  205 DO 210 KORL=KLORLL,K2
      KL=(LORK-1)*ND+KORL
      IF(IW.EQ.2) KL=(KORL-1)*ND+LORK
      KLB=KL+ISHIFT
      Q(KL,1,J)=Q(KLB,1,J)
  210 Q(KL,5,J)=P(KLB,J)/GD
      IF (I1.EQ.1) GO TO 215
      I1=1
      LORK=LORKW+1
      ISHIFT=-ISHIFT
      GO TO 205
  215 IF (JKLW-JLKW) 250,220,230
C     CORNER POINTS AT MUTUAL TRAILING EDGE OF THE TWO WALLS
  220 J=JKLW
      LORK=LORKW
      KORL=KORLW
      I1=0
      ISHIFT=-1-ND
  225 KL=(LORK-1)*ND+KORL
      IF(IW.EQ.2) KL=(KORL-1)*ND+LORK
      KLB=KL+ISHIFT
      Q(KL,1,J)=Q(KLB,1,J)
      Q(KL,5,J)=P(KLB,J)/GD
      IF (I1.EQ.1) GO TO 250
      I1=1
      LORK=LORK+1
      KORL=KORL+1
      ISHIFT=-ISHIFT
      GO TO 225
C     INNER CORNER AND UPPER EDGE OF PROTRUDING SIDEPLATE
C      PROTRUDING EDGE POINTS AT LORKW
  230 J=JLKW
C     INNER CORNER
      LORK=LORKW
      KORL=KORLW
      ISHIFT=-ND
      KL=(LORK-1)*ND+KORL
      IF(IW.EQ.1) GO TO 231
      ISHIFT=-1
      KL=(KORL-1)*ND+LORK
  231 KLB=KL+ISHIFT
      Q(KL,1,J)=Q(KLB,1,J)
      Q(KL,5,J)=P(KLB,J)/GD
C     PROTRUDING EDGE POINTS AT LORKW
      J=JKLW
      LORK=LORKW
      KORL=KORLW
      I1=0
      ISHIFT=-ND
      IF(IW.EQ.2) ISHIFT=-1
  233 KL=(LORK-1)*ND+KORL
      IF(IW.EQ.2) KL=(KORL-1)*ND+LORK
      KLB=KL+ISHIFT
      Q(KL,1,J)=Q(KLB,1,J)
      Q(KL,5,J)=P(KLB,J)/GD
      IF(I1.EQ.1) GO TO 234
      I1=1
      KORL=KORL+1
      GO TO 233
C      UPPER EDGES
  234 LORK=LORKW+1
      KORL=KORLW
      I1=0
  235 KL=(LORK-1)*ND+KORL
      KLB=KL+ND
      IF (IW.EQ.1) GO TO 240
      KL=(KORL-1)*ND+LORK
      KLB=KL+1
  240 DO 245 J=JLKW,JKLW
      Q(KL,1,J)=Q(KLB,1,J)
  245 Q(KL,5,J)=P(KLB,J)/GD
      IF (I1.EQ.1) GO TO 250
      I1=1
      KORL=KORL+1
      GO TO 235
  250 CONTINUE
  255 CONTINUE
C     FILTER SOLUTION AT INFLOW BOUNDARY
      INFILT=0
      IF(RM.EQ.0.) RETURN
      IF(INFILT.EQ.0.OR.J1BC.NE.1) RETURN
      IF(NC.EQ.0) RETURN
      J=1
      L1=1
      IF(LAXIS.EQ.0.AND.L1BC.EQ.1.AND.JL1WL.EQ.1) L1=2
      L2=LM
      IF(LW.GT.0) L2=LW-1
      K1=1
      IF(KPLANE.EQ.0.AND.K1BC.EQ.1.AND.JK1WL.EQ.1) K1=2
      K2=KM
      IF(KW.GT.0) K2=KW-1
      IF(KPLANE.EQ.1) K2=1
C     FILTER DENSITY OVER L
      DO 350 K=K1,K2
      DO 320 L=L1,L2
      KL=(L-1)*ND+K
      KLB=KL-ND
      KLR=KL+ND
      IF(L-1)315,305,315
  305 S(KL,1,J)=.5*(Q(KL,1,J)+Q(KLR,1,J))
      GO TO 320
  315 S(KL,1,J)=.25*(Q(KLB,1,J)+Q(KLR,1,J)+2.*Q(KL,1,J))
  320 CONTINUE
      DO 330 L=L1,L2
  330 Q(KL,1,J)=S(KL,1,J)
  350 CONTINUE
C     FILTER DENSITY OVER K
      DO 500 L=L1,L2
      KL1=(L-1)*ND
      IF(K1.EQ.K2) GO TO 400
      DO 380 K=K1,K2
      KL=KL1+K
      KLB=KL-1
      KLR=KL+1
      IF(K-1) 365,355,365
  355 S(KL,1,J)=.5*(Q(KL,1,J)+Q(KLR,1,J))
      GO TO 380
  365 S(KL,1,J)=.25*(Q(KLB,1,J)+Q(KLR,1,J)+2.*Q(KL,1,J))
  380 CONTINUE
      DO 390 K=K1,K2
      KL=KL1+K
  390 Q(KL,1,J)=S(KL,1,J)
C     RECOMPUTE MOMENTUM FLUXES AND ENERGY IN TERMS OF FILTERED DENSITY
C     AND OF KNOWN TOTAL PRESSURE, TOTAL TEMPERATURE, AND VELOCITY
C     VECTOR INCLINATION
  400 DO 450 K=K1,K2
      KL=KL1+K
      RGMI=(Q(KL,1,J)*HTOT(KL)/PTOT(KL))**GAMI
      WRITE(6,420) K,L,KL,J,Q(KL,1,J),HTOT(KL),PTOT(KL),GAMI
  420 FORMAT(43H K,L,KL,J,Q(KL),1,J),HTOT(KL,PTOT(KL),GAMI=
     1  4I5,1P4E12.5)
      WRITE(6,425) RGMI
  425 FORMAT(7H  RGMI=1PE12.5)
      VSQ=2.*HTOT(KL)*(1.-RGMI)/GAMI
      Q(KL,2,J)=Q(KL,1,J)*SQRT(VSQ/(1.+VOU(KL)**2+WOU(KL)**2))
      Q(KL,3,J)=Q(KL,2,J)*VOU(KL)
      Q(KL,4,J)=Q(KL,2,J)*WOU(KL)
  450 Q(KL,5,J)=Q(KL,1,J)*HTOT(KL)*(GAMMA-GAMI*RGMI)/GD
  500 CONTINUE
      RETURN
      END
*DECK BCALJ
      SUBROUTINE BCALJ(K,L)
*CALL BASE
*CALL RHSBCC
       COMMON/AXISYM/LAXIS
C     SET UP MATRICES FOR ALGEBRAIC BC AT UPPER & LOWER BDYS & INT WALLS
      DO 70 IW=1,2
      CALL SW
      KORL=K
      LORK=L
      IF (IW.EQ.1) GO TO 10
      KORL=L
      LORK=K
   10 IF (LORK.LT.LORKMX) GO TO 35
C     BC AT K=KMAX OR L=LMAX
      IF (LKMXBC.EQ.0.OR.LKMXBC.GE.3) GO TO 35
      IF (LKMXBC-2) 25,15,15
   15 DO 20 J=JB,JMAX
   20 CALL BCFMAT(J,K,L,1)
      GO TO 35
   25 DO 30 J=JB,JMAX
   30 CALL BCWMAT(J,K,L,1)
      IF(JMAXBC.EQ.2) CALL BCJMAT(JMAX,K,L,1)
C     BC AT L=1 OR K=1
   35 IF (LORK.GT.1) GO TO 60
C     TEST FOR INTERMEDIATE WALL NORMAL TO LORK=1 SURFACE
      IF (KORLW) 40,40,55
   40 JBCMIN=MAX0(JLK1WL,JB)
      JBCMAX=JLK1WU
      IF (JBCMAX.LE.JBCMIN) GO TO 60
      DO 50 J=JBCMIN,JBCMAX
   50 CALL BCWMAT(J,K,L,1)
      GO TO 60
   55 IF (KORL.LT.KORLW) GO TO 40
   60 CONTINUE
C     BC AT INTERMEDIATE WALLS
      IF (LORKW.EQ.0) GO TO 70
      IF (LORK.LT.LORKW.OR.LORK.GT.LORKW+1) GO TO 70
      IF (KORLW.GT.0.AND.KORL.GT.KORLW+1) GO TO 70
      JBCMIN=JB
      JBCMAX=JLKW
      IF (JBCMAX.LE.JBCMIN) GO TO 70
      DO 65 J=JBCMIN,JBCMAX
   65 CALL BCWMAT(J,K,L,1)
   70 CONTINUE
C     SET UP MATRIX ELEMENTS FOR USE OF LAGGED B.C.
C     SINGULAR AXIS L= FOR AXIAL SYMMETRY (LAXIS=1)
      IF (L.NE.LAXIS) GO TO 80
      DO 75 J=1,JMAX
   75 CALL ZERODQ(J)
      GO TO 155
C     INTERNAL CORNER AT KMAX,LMAX
   80 IF (KMAXBC.NE.1.OR.LMAXBC.NE.1) GO TO 90
      IF (K*L.NE.KMAX*LMAX) GO TO 90
      DO 85 J=1,JMAX
   85 CALL ZERODQ(J)
      GO TO 155
C     CORNER AT L=1, K=1
   90 IF (K*L.NE.1) GO TO 100
      IF (JK1WU*JL1WU.LE.0) GO TO 100
      J2=MIN0(JK1WU,JL1WU)
      J1=MAX0(1,JK1WL,JL1WL)
      IF (J2.LT.J1) GO TO 100
      DO 95 J=J1,J2
   95 CALL ZERODQ(J)
C     TRAILING EDGES AND CORNERS OF INTERMEDIATE WALLS
      IF(KW.LE.0 .AND. LW.LE.0) GO TO 155
  100 DO 135 IW=1,2
      CALL SW
      KORL=K
      LORK=L
      IF (IW.EQ.1) GO TO 105
      KORL=L
      LORK=K
  105 IF (KORLW.LE.0) GO TO 135
      IF (KORL.LT.KORLW.OR.KORL.GT.KORLW+1) GO TO 135
C     TRAILING EDGE OF WALL KORL=KORLW
      IF (LORKW) 115,115,110
  110 IF (LORK.GT.LORKW) GO TO 125
  115 CALL ZERODQ(JKLW)
C     CORNER OF INTERSECTION WITH PLUG AT LORK=1
      IF (LORK.NE.1.OR.KORL.NE.KORLW) GO TO 125
      J2=JLK1WU
      J1=MAX0(JLK1WL,1)
      IF (J2.LT.J1) GO TO 125
      DO 120 J=J1,J2
  120 CALL ZERODQ(J)
C     INNER EDGE OF PROTRUDING SIDEPLATE
  125 IF(LORKW.EQ.0) GO TO 135
      IF (LORK.NE.LORKW+1.OR.KORL.NE.KORLW) GO TO 135
      IF (JLKW.GE.JKLW) GO TO 135
      DO 130 J=JLKW,JKLW
  130 CALL ZERODQ(J)
  135 CONTINUE
C     CORNERS WHERE THE INTERMEDIATE WALLS INTERSECT
      IF (KW*LW.LE.0) GO TO 155
      IF (L.NE.LW.OR.K.NE.KW) GO TO 145
      J2=MIN0(JKW,JLW)
      DO 140 J=1,J2
  140 CALL ZERODQ(J)
      GO TO 155
  145 IF (L.NE.LW+1.OR.K.NE.KW+1) GO TO 155
      J2=MAX0(JKW,JLW)
      DO 150 J=1,J2
  150 CALL ZERODQ(J)
  155 CONTINUE
C
C     PRESEVE INITIAL DATA AT INFLOW PLANE J=1 IN REGION OUTSIDE NOZZLE
      IF(LW) 2100,2100,2300
 2100 IF(KW.LE.0) GO TO 3000
      IF(K.LE.KW) GO TO 3000
 2200 CALL ZERODQ(1)
      GO TO 3000
 2300 IF(KW) 2400,2400,2500
 2400 IF(L.LE.LW) GO TO 3000
      GO TO 2200
 2500 IF(K.GT.KW .AND. L.GT.LW) CALL ZERODQ(1)
 3000 CONTINUE
C
C     LAGGED B.C. FOR NOZZLE LIP WALL POINTS AT J=JMAX WHEN JMAXBC=2
C       (OUTFLOW B.C. P=FSP IMPOSED)
C
      IF(JMAXBC .NE. 2) GO TO 9999
      IF(ITW .EQ. 1) GO TO 9999
      J = JMAX
C
C     HORIZONTAL WALL AT L=LMAX
C
      IF(L .EQ. LMAX .AND. LMAXBC .EQ. 1) CALL ZERODQ(J)
C
C     VERTICAL WALL AT K=KMAX
C
      IF(KPLANE .EQ. 1) GO TO 9999
      IF(K .EQ. KMAX .AND. KMAXBC .EQ. 1) CALL ZERODQ(J)
 9999 RETURN
      END
*DECK BCALKL
      SUBROUTINE BCALKL(JJ,LORK,IWW)
*CALL BASE
*CALL RHSBCC
      COMMON/AXISYM/LAXIS
      ISW=IWW+1
      IW=IWW
      J=JJ
      CALL SW
C     SET UP MATRICES FOR ALGEBRAIC BC AT LOWER,UPPER BDYS AND INT WALLS
      IF (LORK.LT.LORKMX) GO TO 30
C     BC AT LORK=LORKMX
      IF (LKMXBC.LT.1.OR.LKMXBC.GT.2) GO TO 30
      IF (LKMXBC-2) 20,10,10
C     FREESTREAM BC
   10 DO 15 K=KLORLL,KUORLU
      KR=K
      LR=LORK
      IF (IW.EQ.1) GO TO 15
      KR=LORK
      LR=K
   15 CALL BCFMAT(J,KR,LR,ISW)
      GO TO 30
C     WALL BC
   20 DO 25 K=KLORLL,KUORLU
      KR=K
      LR=LORK
      IF (IW.EQ.1) GO TO 23
      KR=LORK
      LR=K
   23 CALL BCWMAT(J,KR,LR,ISW)
      IF(J.EQ.JMAX .AND. JMAXBC.EQ.2) CALL BCJMAT(J,KR,LR,ISW)
   25 CONTINUE
C     BC AT INTERMEDIATE WALLS
   30 IF (KORLW.LE.0) GO TO 40
      IF (J.GT.JKLW) GO TO 40
      IF (LORKW.GT.0.AND.LORK.GT.LORKW+1) GO TO 50
      KR=KORLW
      LR=LORK
      IF(IW.EQ.1) CALL BCWMAT(J,KR+1,LR,ISW)
      IF (IW.EQ.1) GO TO 35
      KR=LORK
      LR=KORLW
      CALL BCWMAT(J,KR,LR+1,ISW)
   35 CALL BCWMAT(J,KR,LR,ISW)
   40 IF (LORKW.LE.0) GO TO 50
      IF (J.GT.JLKW) GO TO 50
      IF (LORK.LT.LORKW.OR.LORK.GT.LORKW+1) GO TO 50
      KBCMAX=KORLMX
      IF(KORLW.GT.0) KBCMAX=KORLW+1
      DO 45 K=KLORLL,KBCMAX
      KR=K
      LR=LORK
      IF (IW.EQ.1) GO TO 45
      KR=LORK
      LR=K
   45 CALL BCWMAT(J,KR,LR,ISW)
C    BC AT LORK=1
   50 IF (LORK.NE.1.OR.LK1BC.EQ.0) GO TO 60
      KBCMAX=KORLMX
      IF(KORLW.GT.0) KBCMAX=KORLW
      IF (J.LT.JLK1WL.OR.J.GT.JLK1WU) GO TO 65
      DO 55 K=KLORLL,KBCMAX
      KR=K
      LR=LORK
      IF (IW.EQ.1) GO TO 53
      KR=LORK
      LR=K
   53 CALL BCWMAT(J,KR,LR,ISW)
      IF(J.EQ.JMAX .AND. JMAXBC.EQ.2) CALL BCJMAT(J,KR,LR,ISW)
   55 CONTINUE
C    BC AT KORL=1
   60 IF (KL1BC.EQ.0) GO TO 85
C    TEST FOR INTERMEDIATE WALL
   65 IF (LORKW) 70,70,80
   70 IF (J.LT.JKL1WL.OR.J.GT.JKL1WU) GO TO 85
      KR=1
      LR=LORK
      IF (IW.EQ.1) GO TO 75
      KR=LORK
      LR=1
   75 CALL BCWMAT(J,KR,LR,ISW)
      IF(J.EQ.JMAX .AND. JMAXBC.EQ.2) CALL BCJMAT(J,KR,LR,ISW)
      GO TO 85
   80 IF (LORK.LE.LORKW) GO TO 70
C    BC AT KORL=KORLMX
   85 IF (KLMXBC.LT.1.OR.KLMXBC.GT.2) GO TO 105
      KR=KORLMX
      LR=LORK
      IF (IW.EQ.1) GO TO 90
      KR=LORK
      LR=KORLMX
   90 IF (KLMXBC-2) 100,95,95
C    FREESTREAM BC
   95 CALL BCFMAT(J,KR,LR,ISW)
      GO TO 105
C    WALL BC
  100 CALL BCWMAT(J,KR,LR,ISW)
      IF(J.EQ.JMAX .AND. JMAXBC.EQ.2) CALL BCJMAT(J,KR,LR,ISW)
  105 CONTINUE
C     SET UP MATRIX ELEMENTS FOR USE OF LAGGED B.C.
C     SINGULAR AXIS L=1 FOR AXISYMMETRIC FLOW
      IF (IW.NE.1.OR.LORK.NE.LAXIS) GO TO 115
      DO 110 K=KLORLL,KUORLU
  110 CALL ZERODQ(K)
  115 IF(IW.EQ.2.AND.LAXIS.EQ.1) CALL ZERODQ(1)
C
C     NOZZLE LIP WALL POINTS AT J=JMAX FOR JMAXBC=2
C       (OUTFLOW B.C. P=FSP IMPOSED)
C
      IF(JMAXBC .NE. 2 .OR. J .NE. JMAX) GO TO 700
      IF(KPLANE .EQ. 1 .AND. ISW .EQ. 2) GO TO 700
      IF(ITW .EQ. 1) GO TO 700
C
C     WALL AT LORK=LORKMX
C
  630 IF(LORK .NE. LORKMX) GO TO 680
      IF(LKMXBC .NE. 1) GO TO 680
      DO 670 K=KLORLL,KORLMX
  670 CALL ZERODQ(K)
  680 CONTINUE
C
C      WALL AT KORL=KORLMX
C
      IF(KLMXBC .EQ. 1) CALL ZERODQ(KORLMX)
  700 CONTINUE
C     INTERNAL CORNER AT K,L=KMAX,LMAX
      IF (LKMXBC*KLMXBC.NE.1) GO TO 120
      IF (LORK.NE.LORKMX) GO TO 120
      CALL ZERODQ(KORLMX)
C     CORNERS IN PLANE LORK=1
  120 IF (LORK.NE.1) GO TO 125
C     CORNER AT (K,L)=(1,1)
      IF(J.GE.MAX0(JLK1WL,JKL1WL).AND.J.LE.MIN0(JLK1WU,JKL1WU))
     1     CALL ZERODQ(1)
C     CORNER WHERE PLUG INTERSECTS SIDEWALL
      IF (KORLW.LE.0) GO TO 125
      IF(J.GE.JLK1WL.AND.J.LE.JLK1WU) CALL ZERODQ(KORLW)
C     TRAILING EDGES OF LONE INTERMEDIATE WALLS
  125 IF (KORLW.GT.0) GO TO 135
C     TRAILING EDGE OF LONE HORIZONTAL WALL
      IF (LORKW.LE.0.OR.J.NE.JLKW) GO TO 175
      IF (LORK.NE.LORKW.AND.LORK.NE.LORKW+1) GO TO 175
      DO 130 K=KLORLL,KUORLU
  130 CALL ZERODQ(K)
      GO TO 175
C     TRAILING EDGE OF LONE VERTICAL WALL
  135 IF (LORKW.GT.0) GO TO 140
      IF (J.NE.JKLW) GO TO 175
      CALL ZERODQ(KORLW)
      CALL ZERODQ(KORLW+1)
      GO TO 175
C     EDGES AND CORNERS OF COEXISTENT INTERMEDIATE WALLS
  140 IF (LORK.GT.LORKW) GO TO 160
      IF (J.NE.JKLW) GO TO 145
C     TRAILING EDGE OF VERTICAL SIDEPLATE
      CALL ZERODQ(KORLW)
      CALL ZERODQ(KORLW+1)
C     CORNERS AND EDGES ON LOWER SURFACE OF HORIZONTAL WALL
  145 IF (LORK.NE.LORKW) GO TO 175
C     INTERSECTION OF HORIZONTAL WALL AND SIDEPLATE
      IF(J.LE.MIN0(JKLW,JLKW)) CALL ZERODQ(KORLW)
C     INTERSECTION OF HORIZONTAL WALL AND VERTICAL PLUG
      IF(JKL1WL.LE.J.AND.J.LE.JKL1WU) CALL ZERODQ(1)
C     TRAILING EDGE OF HORIZONTAL WALL
      IF (J.NE.JLKW) GO TO 155
      K2=KORLW
      IF(JLKW.GT.JKLW) K2=KORLW+1
      DO 150 K=1,K2
  150 CALL ZERODQ(K)
C     INSIDE OF PROTRUDING SIDEPLATE
  155 IF(JKLW.LE.J.AND.J.LE.JLKW) CALL ZERODQ(KORLW+1)
      GO TO 175
  160 IF (LORK.GT.LORKW+1) GO TO 175
C     CORNERS AND EDGES ON UPPER SURFACE OF HORIZONTAL WALL
C     INTERSECTION WITH VERTICAL WALL
      IF(J.LE.MAX0(JKLW,JLKW))CALL ZERODQ(KORLW+1)
C     TRAILING EDGE
      IF (J.NE.JLKW) GO TO 170
      DO 165 K=1,KORLW
  165 CALL ZERODQ(K)
C     INSIDE EDGE OF PROTRUDING SIDEPLATE
  170 IF(JLKW.LE.J.AND.J.LE.JKLW) CALL ZERODQ(KORLW)
  175 CONTINUE
C     PRESERVE INITIAL DATA AT INFLOW PLANE J=1 IN REGION OUTSIDE NOZZLE
      IF(J.NE.1) GO TO 800
      IF(LORKW) 755,755,775
  755 IF(KORLW.LE.0) GO TO 800
  760 K1=KORLW+1
  765 DO 770 K=K1,KORLMX
  770 CALL ZERODQ(K)
      GO TO 800
  775 IF(LORK-LORKW) 780,780,785
  780 IF(KORLW.LE.0) GO TO 800
      GO TO 760
  785 K1=1
      GO TO 765
  800 CONTINUE
      RETURN
      END
*DECK BCFMAT
      SUBROUTINE BCFMAT(J,K,L,ISW)
*CALL BASE
*CALL VARS
*CALL COUNT
      COMMON/FREESP/FSP,FSRHO
*CALL BTRID
C
C     IMPLICIT MATRIX DEFINITIONS FOR ALGEBRAIC FREESTREAM BC AT LORK=LORKMX
C      WHEN LKMXBC=2
C     ISW=1,2,3 FOR THE ADI SWEEPS IN J,K,L DIRECTIONS,RESP. THE SWEEPS
C   MUST BE PERFORMED IN THE ORDER J,K,L.
C
      KL=(L-1)*ND+K
C
C     ENFORCE FREESTREAM AXIAL VELOCITY U=FSMACH , PRESSURE FSP
C     AND DENSITY FSRHO AT J,K,L.
C
      RR=1.0/Q(KL,1,J)
      U=Q(KL,2,J)*RR
      V=Q(KL,3,J)*RR
      W=Q(KL,4,J)*RR
      XKE=0.5*(U**2+V**2+W**2)
      PJAY=GD*(Q(KL,5,J)-XKE*Q(KL,1,J))
      PFSJAY=FSP
      RFSJAY=FSRHO
      I=L
      IF (ISW-2) 15,20,25
   15 I=J
      GO TO 25
   20 I=K
   25 DO 30 M=1,5
      A(I,1,M)=0.0
      B(I,1,M)=0.0
      C(I,1,M)=0.0
      A(I,2,M)=0.0
      B(I,2,M)=0.0
      C(I,2,M)=0.0
      A(I,5,M)=0.0
      B(I,5,M)=0.0
   30 C(I,5,M)=0.0
      B(I,1,1)=1.0 * Q(KL,6,J)
      B(I,2,1)=-FSMACH * Q(KL,6,J)
      B(I,2,2)=1.0 * Q(KL,6,J)
      B(I,5,1)=XKE * Q(KL,6,J)
      B(I,5,2)=-U * Q(KL,6,J)
      B(I,5,3)=-V * Q(KL,6,J)
      B(I,5,4)=-W * Q(KL,6,J)
      B(I,5,5)=1.0 * Q(KL,6,J)
      IF (ISW-2) 35,40,40
   35 F(I,2)=(FSMACH*Q(KL,1,J)-Q(KL,2,J)) * Q(KL,6,J)
      F(I,5)=((PFSJAY-PJAY)/GD) * Q(KL,6,J)
      F(I,1)=(RFSJAY-Q(KL,1,J)) * Q(KL,6,J)
      RETURN
   40 DO 45 N=1,4
   45 F(I,5)=F(I,5)+B(I,5,N)*F(I,N)/Q(KL,6,J)
      F(I,2)=(B(I,2,1)*F(I,1)+B(I,2,2)*F(I,2))/Q(KL,6,J)
      F(I,1)=(B(I,1,1)*F(I,1))/Q(KL,6,J)
      RETURN
      END
*DECK BCJMAT
      SUBROUTINE BCJMAT(J,K,L,ISW)
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
C     COMMON FREESP CONTAINS FREESTREAM PRESSURE
      COMMON/FREESP/FSP
C     IMPLICIT MATRIX DEFINITIONS FOR ALGEBRAIC B.C. AT J=1,JMAX.
C     ISW=1,2,3 FOR THE ADI SWEEPS IN J,K,L DIRECTIONS,RESP. THE SWEEPS
C   MUST BE PERFORMED IN THE ORDER J,K,L.
      KL=(L-1)*ND+K
      RR=1.0/Q(KL,1,J)
      U=Q(KL,2,J)*RR
      V=Q(KL,3,J)*RR
      W=Q(KL,4,J)*RR
      XKE=0.5*(U**2+V**2+W**2)
      PJAY=GD*(Q(KL,5,J)-XKE*Q(KL,1,J))
      PFSJAY=FSP
      I=L
      IF (ISW-2) 10,15,20
   10 I=J
      GO TO 20
   15 I=K
   20 IF (J-1) 50,50,25
   25 IF (JMAXBC.NE.2) GO TO 145
C   ENFORCE FREESTREAM PRESSURE B.C. P=FSP AT OUTFLOW
C   BOUNDARY J=JMAX (INPUT OPTION JMAXBC=2)
C  THIS ALGEBRAIC B.C. IS USED IN PLACE OF THE ENERGY EQUATION
C  TO MAINTAIN CONSISTENCY WITH THE WALL B.C. AT WALL POINTS.
      IF(L .EQ. LMAX .AND. LMAXBC .EQ. 1) GO TO 26
      IF(K .EQ. KMAX .AND. KMAXBC .EQ. 1) GO TO 26
      GO TO 300
   26 CONTINUE
      DO 30 M=1,5
      A(I,5,M)=0.0
      B(I,5,M)=0.0
   30 C(I,5,M)=0.0
      B(I,5,1)=XKE * Q(KL,6,J)
      B(I,5,2)=-U * Q(KL,6,J)
      B(I,5,3)=-V * Q(KL,6,J)
      B(I,5,4)=-W * Q(KL,6,J)
      B(I,5,5)=1.0 * Q(KL,6,J)
      IF (ISW-2) 35,40,40
   35 F(I,5)=((PFSJAY-PJAY)/GD)*Q(KL,6,J)
      GO TO 145
   40 FSVP=0.0
      DO 45 M=1,5
   45 FSVP=FSVP+B(I,5,M)*F(I,M)
      F(I,5)=FSVP/Q(KL,6,J)
C  ENFORCE WALL TEMPERATURE B.C. AT WALL POINTS OF THE OUTFLOW
C  BOUNDARY WHEN TW.GT.0. THE ALGEBRAIC B.C. RHO*TW=
C  GAMMA*(GAMMA-1.)*EPSLON REPLACES THE CONTINUITY EQUATION
C  IN THIS CASE.
  300 IF(ITW.EQ.0) GO TO 145
      THW=TW/GD
      IF(L-LMAX) 46,47,145
   46 IF(K-KMAX) 145,47,145
   47 DO 200 M=1,5
      A(I,1,M)=0.
      B(I,1,M)=0.
  200 C(I,1,M)=0.
      B(I,1,1)=THW*Q(KL,6,J)
      B(I,1,5)=-Q(KL,6,J)
      F(I,5)=0.
      GO TO 145
C     FOR INPUT OPTION J1BC=1, WE
C   ENFORCE TWO IMPLICIT INFLOW B.C.AT J=1 ON VELOCITY VECTOR DIREC-
C   TION COSINES V/U,W/U. THESE REPLACE THE V,W MOMENTUM EQ'S.,RESP.
C     FOR OPTION J1BC=2, WE COMPUTE V,W FROM THEIR MOMENTUM EQUATIONS.
   50 IF (J1BC.GT.1) GO TO 70
      DO 55 N=3,4
      DO 55 M=1,5
      A(I,N,M)=0.0
      B(I,N,M)=0.0
      C(I,N,M)=0.0
   55 CONTINUE
      B(I,3,2)=-VOU(KL) * Q(KL,6,J)
      B(I,3,3)=1.0 * Q(KL,6,J)
      B(I,4,2)=-WOU(KL) * Q(KL,6,J)
      B(I,4,4)=Q(KL,6,J)
      IF (ISW-2) 60,65,65
   60 F(I,3)=(VOU(KL)*Q(KL,2,J)-Q(KL,3,J))*Q(KL,6,J)
      F(I,4)=(WOU(KL)*Q(KL,2,J)-Q(KL,4,J))*Q(KL,6,J)
      GO TO 70
   65 F(I,3)=(B(I,3,2)*F(I,2)+B(I,3,3)*F(I,3))/Q(KL,6,J)
      F(I,4)=(B(I,4,2)*F(I,2)+B(I,4,4)*F(I,4))/Q(KL,6,J)
C   ENFORCE TWO IMPLICIT INFLOW B.C. ON TOTAL ENTHALPY HTOT AND TOTAL
C   PRESSURE PTOT. J1BC.LT.0 SELECTS AN MOC RELATION TO CLOSE THE
C   SYSTEM. FOR J1BC.GT.0 THE SYSTEM IS CLOSED WITH EITHER THE CONTINU
C   ITY, THE U-MOMENTUM, OR THE ENERGY EQUATION BY DEFINING NX=1,2,OR
C   5, RESP.
C     FOR J1BC=3 THE U-MOMENTUM EQ. IS REPLACED BY THE IMPLICIT RELATION
C     D(U)/D(XI)=0
   70 IF (J1BC) 75,80,80
   75 NX=1
      NPT=2
      NHT=5
C   MOC RELATION FOR B(I,NX,M), C(I,NX,M),F(I,NX) GOES HERE
      GO TO 105
   80 CONTINUE
C   DEFINE NX FOR J1BC.GT.0
      NX=INXBC
      IF (LW.LE.0) GO TO 85
      KR=KMAX
      IF(KW.GT.0) KR=KW
      IF(K.GT.KR) GO TO 85
      IF(L.EQ.LW .OR. L.EQ.LW+1) NX=2
   85 IF (KW.LE.0) GO TO 90
      LR=LMAX
      IF(LW.GT.0) LR=LW
      IF(L.GT.LR) GO TO 90
      IF(K.EQ.KW.OR.K.EQ.KW+1) NX=2
   90 IF(L.EQ.LMAX .AND. LMAXBC.EQ.1) NX=2
      IF(K.EQ.KMAX .AND. KMAXBC.EQ.1) NX=2
      IF(L.EQ.1 .AND. L1BC.EQ.1 .AND. JL1WL.LE.J .AND. J.LE.JL1WU) NX=2
      IF(K.EQ.1 .AND. K1BC.EQ.1 .AND. JK1WL.LE.J .AND. J.LE.JK1WU) NX=2
      NPT=2
      NHT=5
      IF (NX-2) 105,95,100
   95 NPT=1
      GO TO 105
  100 NPT=1
      NHT=2
  105 DO 110 M=1,5
      A(I,NPT,M)=0.0
      B(I,NPT,M)=0.0
      C(I,NPT,M)=0.0
      A(I,NHT,M)=0.0
      B(I,NHT,M)=0.0
  110 C(I,NHT,M)=0.0
      KL=KL
      HG=HTOT(KL)/GAMI
      IF (J1BC.EQ.3) GO TO 115
      GCR=HG*(Q(KL,1,J)*HTOT(KL)/PTOT(KL))**GAMI
      B(I,NPT,1)=(XKE-GCR)*Q(KL,6,J)
      B(I,NPT,2)=-U * Q(KL,6,J)
      B(I,NPT,3)=-V * Q(KL,6,J)
      B(I,NPT,4)=-W * Q(KL,6,J)
      B(I,NPT,5)=1.0 * Q(KL,6,J)
      GO TO 125
  115 IF (ISW.GE.2) GO TO 120
      B(I,NPT,1)=Q(KL,2,J)/Q(KL,1,J)**2
      B(I,NPT,2)=-1.0/Q(KL,1,J)
      C(I,NPT,1)=-Q(KL,2,J+1)/Q(KL,1,J+1)**2
      C(I,NPT,2)=1.0/Q(KL,1,J+1)
      GO TO 125
  120 B(I,NPT,2)=1.0
  125 B(I,NHT,1)=(HG-GAMI*XKE)*Q(KL,6,J)
      B(I,NHT,2)=GAMI*U * Q(KL,6,J)
      B(I,NHT,3)=GAMI*V * Q(KL,6,J)
      B(I,NHT,4)=GAMI*W * Q(KL,6,J)
      B(I,NHT,5)=-GAMMA * Q(KL,6,J)
      IF (ISW-2) 130,135,135
  130 IF(J1BC.LE.2) F(I,NPT)=((XKE+(GCR/GAMMA))*Q(KL,1,J)-Q(KL,5,J))
     1  *Q(KL,6,J)
      IF(J1BC.EQ.3)F(I,NPT)=(Q(KL,2,J)/Q(KL,1,J)-Q(KL,2,J+1)/Q(KL,1,J+1)
     1  )*Q(KL,6,J)
      F(I,NHT)=(GAMMA*Q(KL,5,J)-Q(KL,1,J)*(HG+GAMI*XKE))*Q(KL,6,J)
      GO TO 145
  135 FSVP=0.0
      FSVH=0.0
      DO 140 M=1,5
      FSVP=FSVP+B(I,NPT,M)*F(I,M)
  140 FSVH=FSVH+B(I,NHT,M)*F(I,M)
      F(I,NPT)=FSVP/Q(KL,6,J)
      F(I,NHT)=FSVH/Q(KL,6,J)
  145 RETURN
      END
*DECK BCLK1
      SUBROUTINE BCLK1(IWW,JJ,K)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL RHSBCC
C   IMPLICIT WALL BC AT GRID POINTS LK=1 AND J.GE.JLK1WL.OR.J.LE.JLK1WU
C     AND SYMMETRY BC OUTSIDE THAT J INTERVAL INPLICIT WALL BC AT
C     INTERMEDIATE WALL UPPER SURFACE L=LORKW+1 AND J.LE.JLKW.
C   FILLS IN FIRST ROW OF COEFFICIENT MATRIX TO ACCOUNT FOR IMPLICIT
C   PORTION OF INVISCID FLUX VECTOR TERMS, AND FILLS IN FIRST ROW OF
C   RIGHT HAND SIDE VECTOR AS COMPUTED IN SUBROUTINE RHS. FOR
C   ADIABATIC WALL, THE IMPLICIT PORTION OF THE VISCOUS TERMS IN THE
C   ENERGY EQ FOR ZETA(L) OR ETA(K) DIRECTION ARE ADDED BY SUBR VISMAT.
C     THE EXTRA MATRIX MULTIPLY THAT ACCOUNTS FOR ZERO VELOCITY B.C.
C   AND FOR PRESCRIBED WALL TEMP BC (ITW=1) IS PERFORMED IN  SUBR BCWMAT
      IW=IWW
      J=JJ
      CALL SW
      LKWR=LORKW+1
      KLWR=KORLW+1
      INTWAL=0
      IF(LORKW.GT.0.AND.J.LE.JLKW) INTWAL=1
      IF(KORLW.GT.0.AND.K.GE.KLWR) INTWAL=0
C     BC AT L=1
      L=1
   50 KL=(L-1)*ND+K
      KLF=KL+ND
      IF(IW.EQ.1) GO TO 75
      KL=(K-1)*ND+L
      KLF=KL+1
   75 RJ=Q(KL,6,J)
      RF=RM*(RJ+Q(KLF,6,J))
C     BC AT L=1
      IF(L.GT.1) GO TO 150
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1 SURFACE
      IF(KORLW) 100,100,130
C     TEST FOR WALL AT L=1
  100 IF(JLK1WL.LE.J.AND.J.LE.JLK1WU) GO TO 150
C     SYMMETRY
      GO TO 190
C     TEST FOR WALL AT L=1
  130 IF(K-KORLW) 100,100,190
  150 DO 180 N=1,5
      DO 170 M=1,5
      A(L,N,M)=0.0
      B(L,N,M)=-2.0*D(L,N,M)
  170 C(L,N,M)=2.0*D(L+1,N,M)
      B(L,N,N)=B(L,N,N)+(WTFAC*RJ+RF)
      C(L,N,N)=C(L,N,N)+(-RF+(1.0-WTFAC)*RJ)
  180 F(L,N)=S(KL,N,J)
      IF(L.EQ.1) GO TO 250
C
C   SYM BC J.LT.JLK1WL.OR.J.GT.JLK1WU FOR INVISCID TERMS AT LWER BDY L=1
C
  190 IF(L.GT.1) GO TO 300
      DO 210 N=1,5
      DO 200 M=1,5
      A(L,N,M)=0.0
      B(L,N,M)=0.0
  200 C(L,N,M)=2.0*D(L+1,N,M)
      B(L,N,N)=Q(KL,6,J)+RF
      C(L,N,N)=C(L,N,N)-RF
  210 F(L,N)=S(KL,N,J)
      IF(IW.EQ.2) GO TO 230
      DO 220 M=1,5
  220 C(L,4,M)=0.0
      F(L,4)=0.0
      GO TO 250
  230 DO 240 M=1,5
  240 C(L,3,M)=0.0
      F(L,3)=0.0
C     BC AT INTERMEDIATE WALL
  250 IF(INTWAL.EQ.0) GO TO 300
      L=LKWR
      GO TO 50
  300 RETURN
      END
*DECK BCLKMX
      SUBROUTINE BCLKMX(IWW,JJ,K)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL RHSBCC
C     IMPLICIT OUTFLOW BOUNDARY CONDITIONS (LORKMXBC=2 OR 5),OR WALL BC
C   (LKMXBC=1), OR SYMMETRY BC (LKMXBC=3,4) AT GRID POINTS L=LORKMX.
C     IMPLICIT WALL BC INTERMEDIATE WALL LOWER SURFACE L=LORKW,J.LE.JLKW
C     FILLS IN LAST ROW OF BLOCK - TRIDIAGONAL COEFFICIENT MATRIX TO
C    ACCOUNT FOR IMPLICIT PORTION OF INVISCID FLUX VECTOR TERMS, AND
C   FILLS INLAST ROW OF RIGHT HAND SIDE VECTOR AS COMPUTED IN SUBR.
C   RHS. THE IMPLICIT PORTION OF THE VISCOUS TERMS FOR THE ZETA(L)
C    OR ETA(K)  DIRECTION ARE ADDED IN BY SUBROUTINE VISMAT.
      IW=IWW
      J=JJ
      CALL SW
      LKWR=LORKW+1
      KLWR=KORLW+1
      INTWAL=0
      IF(LORKW.GT.0.AND.J.LE.JLKW) INTWAL=1
      IF(KORLW.GT.0.AND.K.GE.KLWR) INTWAL=0
C     BC AT UPPER BOUNDARY LORKMX
      L=LORKMX
  100 KL=(L-1)*ND+K
      KLR=KL-ND
      IF(IW.EQ.1) GO TO 200
      KL=(K-1)*ND+L
      KLR=KL-1
  200 CONTINUE
      RJ=Q(KL,6,J)
      RR = RM * (RJ+Q(KLR,6,J))
      IF(L.EQ.LORKMX.AND.(LKMXBC.EQ.3.OR.LKMXBC.EQ.4)) GO TO 500
      DO 400 N=1,5
      DO 300 M=1,5
      C(L,N,M)=0.0
      A(L,N,M) = -2.0*D(L-1,N,M)
  300 B(L,N,M) = 2.0*D(L,N,M)
      A(L,N,N) = A(L,N,N) +(-RR+(1.0-WTFAC)*RJ)
      B(L,N,N) = WTFAC*RJ + RR + B(L,N,N)
  400 F(L,N) = S(KL,N,J)
C     SYMMETRY BC FOR INVISCID TERMS AT UPPER BOUNDARY LORKMX
      IF(L-LORKMX) 900,850,850
  500 DO 700 N=1,5
      DO 600 M=1,5
      C(L,N,M)=0.0
      B(L,N,M)=0.0
  600 A(L,N,M)=-2.0*D(L-1,N,M)
      B(L,N,N)=Q(KL,6,J)+RR
      A(L,N,N)=A(L,N,N)-RR
  700 F(L,N)=S(KL,N,J)
      DO 800 M=1,5
      A(L,LKMXBC,M)=0.0
  800 F(L,LKMXBC)=0.0
C     BC AT INTERMEDIATE WALL
  850 IF(INTWAL.EQ.0) GO TO 900
      L=LORKW
      GO TO 100
  900 RETURN
      END
*DECK BCWMAT
      SUBROUTINE BCWMAT(J,K,L,ISW)
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
      MDOT=1
C
C     IMPLICIT MATRIX DEFINITIONS FOR ALGEBRAIC WALL BC AT (J,K,L)
C     ISW=1,2,3 FOR THE ADI SWEEPS IN J,K,L DIRECTIONS,RESP. THE SWEEPS
C   MUST BE PERFORMED IN THE ORDER J,K,L.
C
      KL=(L-1)*ND+K
      IF(ITW) 10,10,40
C     VELOCITY B.C. FOR ADIABATIC WALL AND ZERO MASS TRANSFER
   10 I=L
      IF (ISW-2) 15,20,25
   15 I=J
      GO TO 25
   20 I=K
   25 CONTINUE
      DO 30 N=2,4
      DO 30 M=1,5
      A(I,N,M)=0.0
      B(I,N,M)=0.0
      C(I,N,M)=0.0
   30 CONTINUE
      DO 35 N=2,4
      F(I,N)=0.0
   35 B(I,N,N)=Q(KL,6,J)
      RETURN
C?VELOCITY AND TEMPERATURE B.C. FOR PRESCRIBED WALL TEMP. TW AND
C?NORMAL MASS INJECTION RATE MWDOT
C?ZETA(L) DIRECTION ISASSUMED NORMAL TO WALL
C      WMDOT IS ASSUMED ZERO BELOW
C
   40 THW=TW/GD
      I=L
      IF(ISW-2) 45,50,55
   45 I=J
      GO TO 55
   50 I=K
  55  DO 65 N=2,5
      DO 60 M=1,5
      A(I,N,M)=0.
      B(I,N,M)=0.
   60 C(I,N,M)=0.
   65 B(I,N,N)=Q(KL,6,J)
      B(I,5,1)=THW*Q(KL,6,J)
      B(I,5,5)=-Q(KL,6,J)
      DO 70 N=2,5
   70 F(I,N)=0.0
      RETURN
      END
*DECK BCXIN
      SUBROUTINE BCXIN(KL)
*CALL BASE
*CALL VARS
*CALL BTRID
C    IMPLICIT INFLOW B.C. AT GRID POINTS J=1
C   FILLS IN FIRST ROW OF COEFFICIENT MATRIX TO ACCOUNT FOR IMPLICIT
C   PORTION OF INVISCID FLUX VECTOR TERMS, AND FILLS IN FIRST ROW OF
C   RIGHT HAND SIDE VECTOR AS COMPUTED IN SUBRS. RHS AND VISRHS.
C     THE EXTRA MATRIX MULTIPLY THAT ACCOUNTS FOR ALGEBRAIC B.C.
C   IS PERFORMED IN SUBR. BCJMAT
C
      J=1
      JF=J+1
      RJ=Q(KL,6,J)
      RF=RM *(RJ+Q(KL,6,JF))
      DO 20 N=1,5
      DO 15 M=1,5
      A(J,N,M)=0.0
      B(J,N,M)=-2.0*D(J,N,M)
   15 C(J,N,M)=2.0*D(J+1,N,M)
      B(J,N,N)=B(J,N,N)+(WTFAC*RJ+RF)
      C(J,N,N)=C(J,N,N)+(-RF+(1.0-WTFAC)*RJ)
   20 F(J,N)=S(KL,N,J)
   25 RETURN
      END
*DECK BCXOUT
      SUBROUTINE BCXOUT(KL)
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
C   IMPLICIT OUTFLOW BOUNDARY CONDITIONS IN XI DIRECTION
C   FILL IN LAST ROW OF BLOCK TRIDIAGONAL COEFF MATRIX WITH
C   EQUATION BETWEEN DELTAQ(JMAX) AND DELTAQ(JM).
C   EQUATION IS FIRST ORDER BACKWARD DIFFERENCING IN XI OF
C   THE P.D.E. AT J=JMAX.
C   RIGHT HAND SIDE OF EQUATION IS COMPUTED IN SUBR. RHS
      J=JMAX
      RJ=Q(KL,6,J)
      RR=RM * (RJ+Q(KL,6,J-1))
      DO 20 N=1,5
      DO 15 M=1,5
      A(J,N,M)=-2.0*D(J-1,N,M)
   15 B(J,N,M)=2.0*D(J,N,M)
      A(J,N,N)=A(J,N,N)-RR + (1.0-WTFAC)*RJ
      B(J,N,N)=(WTFAC*RJ+RR)+B(J,N,N)
   20 F(J,N)=S(KL,N,J)
      RETURN
      END
*DECK BTRI
      SUBROUTINE BTRI(IL,IU)
*CALL BTRID
      DIMENSION H(5,5)
      REAL L11,L21,L22,L31,L32,L33,L41,L42,L43,L44,L51,L52,L53,L54,L55
      IS=IL+1
      IE=IU-1
C     INSERT LUDEC
      L11=1./B(IL,1,1)
      L21=B(IL,2,1)
      U12=B(IL,1,2)*L11
      L22=1./(B(IL,2,2)-L21*U12)
      U13=B(IL,1,3)*L11
      U14 = B(IL,1,4)*L11
      U15=B(IL,1,5)*L11
      L31=B(IL,3,1)
      L32=B(IL,3,2)-L31*U12
      U23=(B(IL,2,3)-L21*U13)*L22
      L33=1./(B(IL,3,3)-U13*L31-U23*L32)
      U24=(B(IL,2,4)-L21*U14)*L22
      U25=(B(IL,2,5)-L21*U15)*L22
      L41=B(IL,4,1)
      L42=B(IL,4,2)-L41*U12
      L43=B(IL,4,3)-L41*U13-L42*U23
      U34=(B(IL,3,4)-L31*U14-L32*U24)*L33
      L44=1./(B(IL,4,4)-U14*L41-U24*L42-U34*L43)
      U35=(B(IL,3,5)-L31*U15-L32*U25)*L33
      L51=B(IL,5,1)
      L52=B(IL,5,2)-L51*U12
      L53=B(IL,5,3)-L51*U13-L52*U23
      L54=B(IL,5,4)-L51*U14-L52*U24-L53*U34
      U45=(B(IL,4,5)-L41*U15-L42*U25-L43*U35)*L44
      L55=1./(B(IL,5,5)-L51*U15-L52*U25-L53*U35-L54*U45)
C     COMPUTE LITTLE R S
      D1=L11*F(IL,1)
      D2=L22*(F(IL,2)-L21*D1)
      D3=L33*(F(IL,3)-L31*D1-L32*D2)
      D4=L44*(F(IL,4)-L41*D1-L42*D2-L43*D3)
      D5=L55*(F(IL,5)-L51*D1-L52*D2-L53*D3-L54*D4)
C     COMPUTE BIG R S
      F(IL,5)=D5
      F(IL,4)=D4-U45*D5
      F(IL,3)=D3-U34*F(IL,4)-U35*D5
      F(IL,2)=D2-U23*F(IL,3)-U24*F(IL,4)-U25*D5
      F(IL,1)=D1-U12*F(IL,2)-U13*F(IL,3)-U14*F(IL,4)-U15*D5
C     COMPUTE C PRIME FOR FIRST ROW
      DO 10 M=1,5
      D1=L11*C(IL,1,M)
      D2=L22*(C(IL,2,M)-L21*D1)
      D3=L33*(C(IL,3,M)-L31*D1-L32*D2)
      D4=L44*(C(IL,4,M)-L41*D1-L42*D2-L43*D3)
      D5=L55*(C(IL,5,M)-L51*D1-L52*D2-L53*D3-L54*D4)
      B(IL,5,M)=D5
      B(IL,4,M)=D4-U45*D5
      B(IL,3,M) = D3-U34*B(IL,4,M)-U35*D5
      B(IL,2,M) = D2-U23*B(IL,3,M)-U24*B(IL,4,M)-U25*D5
   10 B(IL,1,M) = D1-U12*B(IL,2,M)-U13*B(IL,3,M)-U14*B(IL,4,M)-U15*D5
      DO 30 I=IS,IE
      IR=I-1
C      COMPUTE B PRIME*BIGR
      DO 15 N=1,5
   15 F(I,N)=F(I,N)-A(I,N,1)*F(IR,1)-A(I,N,2)*F(IR,2)-A(I,N,3)*F(IR,3)
     1 -A(I,N,4)*F(IR,4)-A(I,N,5)*F(IR,5)
C      COMPUTE B PRIME
      DO 20 M=1,5
      DO 20 N=1,5
   20 H(N,M)=B(I,N,M)-A(I,N,1)*B(IR,1,M)-A(I,N,2)*B(IR,2,M)-A(I,N,3)
     1 *B(IR,3,M)-A(I,N,4)*B(IR,4,M)-A(I,N,5)*B(IR,5,M)
C     INSERT LUDEC AGAIN
      L11=1./H(1,1)
      L21=H(2,1)
      U12=H(1,2)*L11
      L22=1./(H(2,2)-L21*U12)
      U13=H(1,3)*L11
      U14=H(1,4)*L11
      U15=H(1,5)*L11
      L31=H(3,1)
      L32=H(3,2)-L31*U12
      U23=(H(2,3)-L21*U13)*L22
      L33=1./(H(3,3)-U13*L31-U23*L32)
      U24=(H(2,4)-L21*U14)*L22
      U25=(H(2,5)-L21*U15)*L22
      L41=H(4,1)
      L42=H(4,2)-L41*U12
      L43=H(4,3)-L41*U13-L42*U23
      U34=(H(3,4)-L31*U14-L32*U24)*L33
      L44=1./(H(4,4)-U14*L41-U24*L42-U34*L43)
      U35=(H(3,5)-L31*U15-L32*U25)*L33
      L51=H(5,1)
      L52=H(5,2)-L51*U12
      L53=H(5,3)-L51*U13-L52*U23
      L54=H(5,4)-L51*U14-L52*U24-L53*U34
      U45=(H(4,5)-L41*U15-L42*U25-L43*U35)*L44
      L55=1./(H(5,5)-L51*U15-L52*U25-L53*U35-L54*U45)
C     COMPUTE LITTLE R#S
      D1=L11*F(I,1)
      D2=L22*(F(I,2)-L21*D1)
      D3=L33*(F(I,3)-L31*D1-L32*D2)
      D4=L44*(F(I,4)-L41*D1-L42*D2-L43*D3)
      D5=L55*(F(I,5)-L51*D1-L52*D2-L53*D3-L54*D4)
C     COMPUTE BIG R#S
      F(I,5)=D5
      F(I,4)=D4-U45*D5
      F(I,3)=D3-U34*F(I,4)-U35*D5
      F(I,2)=D2-U23*F(I,3)-U24*F(I,4)-U25*D5
      F(I,1)=D1-U12*F(I,2)-U13*F(I,3)-U14*F(I,4)-U15*D5
C     COMPUTE C PRIMES
      DO 25 M=1,5
      D1=L11*C(I,1,M)
      D2=L22*(C(I,2,M)-L21*D1)
      D3=L33*(C(I,3,M)-L31*D1-L32*D2)
      D4=L44*(C(I,4,M)-L41*D1-L42*D2-L43*D3)
      D5=L55*(C(I,5,M)-L51*D1-L52*D2-L53*D3-L54*D4)
      B(I,5,M)=D5
      B(I,4,M)=D4-U45*D5
      B(I,3,M) = D3-U34*B(I,4,M)-U35*D5
      B(I,2,M) = D2-U23*B(I,3,M)-U24*B(I,4,M)-U25*D5
   25 B(I,1,M) = D1-U12*B(I,2,M)-U13*B(I,3,M)-U14*B(I,4,M)-U15*D5
   30 CONTINUE
      I=IU
      IR=I-1
C      COMPUTE B PRIME*BIG R FOR LAST ROW
      DO 35 N=1,5
   35 F(I,N)=F(I,N)-A(I,N,1)*F(IR,1)-A(I,N,2)*F(IR,2)-A(I,N,3)*F(IR,3)
     1 -A(I,N,4)*F(IR,4)-A(I,N,5)*F(IR,5)
C      COMPUTE B PRIME
      DO 40 M=1,5
      DO 40 N=1,5
   40 H(N,M)=B(I,N,M)-A(I,N,1)*B(IR,1,M)-A(I,N,2)*B(IR,2,M)-A(I,N,3)
     1 *B(IR,3,M)-A(I,N,4)*B(IR,4,M)-A(I,N,5)*B(IR,5,M)
      L11=1./H(1,1)
      L21=H(2,1)
      U12=H(1,2)*L11
      L22=1./(H(2,2)-L21*U12)
      U13=H(1,3)*L11
      U14=H(1,4)*L11
      U15=H(1,5)*L11
      L31=H(3,1)
      L32=H(3,2)-L31*U12
      U23=(H(2,3)-L21*U13)*L22
      L33=1./(H(3,3)-U13*L31-U23*L32)
      U24=(H(2,4)-L21*U14)*L22
      U25=(H(2,5)-L21*U15)*L22
      L41=H(4,1)
      L42=H(4,2)-L41*U12
      L43=H(4,3)-L41*U13-L42*U23
      U34=(H(3,4)-L31*U14-L32*U24)*L33
      L44=1./(H(4,4)-U14*L41-U24*L42-U34*L43)
      U35=(H(3,5)-L31*U15-L32*U25)*L33
      L51=H(5,1)
      L52=H(5,2)-L51*U12
      L53=H(5,3)-L51*U13-L52*U23
      L54=H(5,4)-L51*U14-L52*U24-L53*U34
      U45=(H(4,5)-L41*U15-L42*U25-L43*U35)*L44
      L55=1./(H(5,5)-L51*U15-L52*U25-L53*U35-L54*U45)
      D1=L11*F(I,1)
      D2=L22*(F(I,2)-L21*D1)
      D3=L33*(F(I,3)-L31*D1-L32*D2)
      D4=L44*(F(I,4)-L41*D1-L42*D2-L43*D3)
      D5=L55*(F(I,5)-L51*D1-L52*D2-L53*D3-L54*D4)
      F(I,5)=D5
      F(I,4)=D4-U45*D5
      F(I,3)=D3-U34*F(I,4)-U35*D5
      F(I,2)=D2-U23*F(I,3)-U24*F(I,4)-U25*D5
      F(I,1)=D1-U12*F(I,2)-U13*F(I,3)-U14*F(I,4)-U15*D5
      I=IU
   45 I=I-1
      IP=I+1
      DO 50 N=1,5
   50 F(I,N)=F(I,N)-F(IP,1)*B(I,N,1)-F(IP,2)*B(I,N,2)-F(IP,3)*B(I,N,3)-
     1 F(IP,4)*B(I,N,4)-F(IP,5)*B(I,N,5)
      IF (I.GT.IL) GO TO 45
      RETURN
      END
*DECK DIFFER
      SUBROUTINE DIFFER(RD,IMIN,IMAX,ISYMIN,ISYMAX)
*CALL BTRID
      I1=IMIN+1
      I2=IMAX-1
      DO 100 N=1,5
      DO 100 I=I1,I2
  100 F(I,N)=RD*(C(I+1,N,1)-C(I-1,N,1))
      IF(ISYMIN.GT.0) GO TO 300
      DO 200 N=1,5
  200 F(IMIN,N)=2.0*RD*(C(I1,N,1)-C(IMIN,N,1))
      GO TO 500
  300 DO 400 N=1,5
  400 F(IMIN,N)=2.0*RD*C(I1,N,1)
      F(IMIN,ISYMIN)=0.0
  500 IF(ISYMAX.GT.0) GO TO 700
      DO 600 N=1,5
  600 F(IMAX,N)=2.0*RD*(C(IMAX,N,1)-C(I2,N,1))
      RETURN
  700 DO 800 N=1,5
  800 F(IMAX,N)=-2.0*RD*C(I2,N,1)
      IF(ISYMAX.LT.6) F(IMAX,ISYMAX)=0.0
      RETURN
      END
*DECK DJMET
      SUBROUTINE DJMET(J,KL,XJ,YJ,ZJ)
*CALL BASE
*CALL VARS
      DX2 = .5/DX1
      FAC=1.0
      JP = J+1
      JR = J-1
      IF(J.NE.1) GO TO 100
C     FORWARD DIFFERENCE
      JR=J
      FAC=2.0
      GO TO 200
  100 IF(J.NE.JMAX) GO TO 200
C     BACKWARD DIFFERENCE
      JP=J
      FAC=2.0
  200 XJ = (X(KL,JP)-X(KL,JR))*DX2*FAC
      YJ = (Y(KL,JP)-Y(KL,JR))*DX2*FAC
      ZJ = (Z(KL,JP)-Z(KL,JR))*DX2*FAC
      RETURN
      END
*DECK DKMET
      SUBROUTINE DKMET(J,K,L,KL,XK,YK,ZK)
*CALL BASE
*CALL VARS
C     INTERMEDIATE WALL IN K DIRECTION
      DY2 = .5/DY1
      KP = KL+1
      KR = KL-1
      IF(KW.LE.0.OR.K.LT.KW.OR.K.GT.KW+1.OR.J.GT.JKW) GO TO 200
      IF(LW.GT.0.AND.L.GT.LW) GO TO 200
      IF(K.EQ.KW) GO TO 700
      IF(K.EQ.KW+1) GO TO 500
  200 IF(K.NE.1) GO TO 600
C     TEST FOR PLANAR SYMMETRY
      IF(KPLANE.EQ.0) GO TO 250
      IF(KMAX.EQ.1) GO TO 1000
C     TEST FOR INTERMEDIATE WALL NORMAL TO K=1 SURFACE
  250 IF(LW) 300,300,450
  300 IF(JK1WL.LE.J.AND.J.LE.JK1WU) GO TO 500
C     SYMMETRY
  400 XK=0.0
      ZK=0.0
      YK=2.0*(Y(KP,J)-Y(KL,J))*DY2
      RETURN
C     TEST FOR WALL NORMAL TO K=1 SURFACE
  450 IF(L-LW) 300,300,400
C     FOWARD DIFFERECNE
  500 FAC=2.0
      KR=KL
      GO TO 900
  600 IF(K.NE.KMAX) GO TO 800
C     TEST FOR PLANAR SYMMETRY
      IF(KPLANE.EQ.0) GO TO 650
      IF(KMAX.EQ.1) GO TO 1000
  650 IF(KMAXBC.LT.3.OR.KMAXBC.GT.4) GO TO 700
C     SYMMETRY
      XK=0.0
      YK=0.0
      ZK=0.0
      KR=KR
      IF(KMAXBC.EQ.3) YK=2.0*(Y(KL,J)-Y(KR,J))*DY2
      IF(KMAXBC.EQ.4) ZK=2.0*(Z(KL,J)-Z(KR,J))*DY2
      RETURN
C     BACKWARD DIFFERENCE
  700 KP=KL
      FAC=2.0
      GO TO 900
C     CENTRAL DIFFERENCE
  800 FAC=1.0
  900 KR=KR
      KP=KP
      XK=(X(KP,J)-X(KR,J))*DY2*FAC
      YK=(Y(KP,J)-Y(KR,J))*DY2*FAC
      ZK=(Z(KP,J)-Z(KR,J))*DY2*FAC
      RETURN
C     PLANAR SYMMETRY
 1000 XK=0.0
      YK=1.0
      ZK=0.0
      RETURN
      END
*DECK DLMET
      SUBROUTINE DLMET(J,K,L,KL,XL,YL,ZL)
*CALL BASE
*CALL VARS
      COMMON/AXISYM/LAXIS
C     INTERMEDIATE WALL IN L DIRECTION
      DZ2 = .5/DZ1
      LP = KL+ND
      LR = KL-ND
      IF(LW.LE.0.OR.L.LT.LW.OR.L.GT.LW+1.OR.J.GT.JLW) GO TO 100
      IF(KW.GT.0.AND.K.GT.KW) GO TO 100
      IF(L.EQ.LW) GO TO 700
      IF(L.EQ.LW+1) GO TO 500
  100 IF(L.NE.1) GO TO 600
C     AXIS OF SYMMETRY
      IF(LAXIS.NE.1) GO TO 150
      XL=0.0
      YL=Y(LP,J)-Y(KL,J)
      ZL=Z(LP,J)-Z(KL,J)
      GO TO 999
  150 CONTINUE
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1 SURFACE
      IF(KW) 200,200,400
C     TEST FOR WALL AT L=1
  200 IF(JL1WL.LE.J.AND.J.LE.JL1WU) GO TO 500
C     SYMMETRY
  300 XL=0.0
      YL=0.0
      ZL=2.0*(Z(LP,J)-Z(KL,J))*DZ2
      GO TO 999
C     TEST FOR WALL AT L=1
  400 IF(K-KW) 200,200,300
C     FORWARD DIFFERENCE
  500 LR=KL
      FAC=2.0
      GO TO 900
  600 IF(L.NE.LMAX) GO TO 800
      IF(LMAXBC.LT.3.OR.LMAXBC.GT.4) GO TO 700
C     SYMMETRY
      XL=0.0
      YL=0.0
      ZL=0.0
      LR=LR
      IF(LMAXBC.EQ.4) ZL=2.0*(Z(KL,J)-Z(LR,J))*DZ2
      IF(LMAXBC.EQ.3) YL=2.0*(Y(KL,J)-Y(LR,J))*DZ2
      GO TO 999
C     BACKWARD DIFFERENCE
  700 LP=KL
      FAC=2.0
      GO TO 900
C     CENTRAL DIFFERENCE
  800 FAC=1.0
  900 LR=LR
      LP=LP
      XL=(X(LP,J)-X(LR,J))*DZ2*FAC
      YL=(Y(LP,J)-Y(LR,J))*DZ2*FAC
      ZL=(Z(LP,J)-Z(LR,J))*DZ2*FAC
  999 RETURN
      END
*DECK DYZ
      SUBROUTINE DYZ
*CALL BTRID
*CALL RHSBCC
      I1=I+1
      IF(I.EQ.3) I1=1
      DO 100 L=1,LORKMX
      DO 100 N=1,5
  100 F(L,N)=C(L,N,I)-C(L,N,I1)
      RETURN
      END
*DECK DYZBC1
      SUBROUTINE DYZBC1(ISW)
*CALL BASE
*CALL RHSBCC
*CALL BTRID
C     EVALUATE VISCOUS CROSS TERMS AT A LOWER BDRY K=1 (L=1) WHEN ISW=1,
C     OR AT THE OUTER SURFACE K=KW+1 (L=LW+1) OF AN INTERNAL WALL
C     WHEN ISW=2
      L1=1
      IF(ISW.EQ.2) GO TO 100
      IF(KL1BC.EQ.0) RETURN
C     BOUNDARY AT K=1(L=1)
      IF(J.LT.JKL1WL.OR.J.GT.JKL1WU) GO TO 300
  100 L2=LORKMX
      IF(LORKW.GT.0) L2=LORKW+ISW-1
C     INTERNAL WALL
      IF(ISW.EQ.1) GO TO 150
      IF(J-JKLW) 150,150,110
  110 IF(L2.EQ.LORKMX) RETURN
      IF(J.GT.JLKW) RETURN
      L1=L2-1
C     ADIABATIC WALL BC ON ENERGY EQUATION
  150 I1=I-1
      IF(I.EQ.1) I1=3
      DO 200 L=1,L2
      DO 190 N=1,4
  190 F(L,N)=0.
  200 F(L,5)=2.0*(C(L,5,I)-C(L,5,I1))
      IF(ISW.EQ.2) RETURN
  250 IF(L2.EQ.LORKMX) RETURN
      L1=L2+1
C     SYM BC AT KL=1.AND.J.LT.JKL1WL.OR.J.GT.JKL1WU.AND.LK.GT.LORKW
  300 DO 500 L=L1,LORKMX
      DO 400 N=1,5
  400 F(L,N)=2.0*(C(L,N,2)-C(L,N,1))
  500 F(L,IW+2)=0.0
      RETURN
      END
*DECK DYZBCX
      SUBROUTINE DYZBCX(ISW)
*CALL BASE
*CALL RHSBCC
*CALL BTRID
C     EVALUATE VISCOUS CROSS TERMS AT UPPER BOUNDARY KMAX(LMAX WHEN ISW=1,
C     OR AT THE INNER SURFACE K=KW (L=LW) OF AN INTERNAL WALL WHEN ISW=2
      I1=I-1
      IF(I.EQ.1) I1=3
      L2=LORKMX
      IF(ISW.EQ.2) GO TO 300
      IF(KLMXBC.EQ.0) RETURN
C     BOUNDARY AT K=KMAX(L=LMAX)
C     WALL B.C.(LKMXBC=1),
C     FREE STREAM BC(KLMXBC=2) OR OUTFLOW BC(KLMXBC=5)
   80 DO 90 L=1,L2
   90 F(L,1)=0.0
      DO 100 N=2,5
      DO 100 L=1,L2
  100 F(L,N)=2.0*(C(L,N,I)-C(L,N,I1))
      IF(ISW.EQ.2) RETURN
      IF(KLMXBC.EQ.3.OR.KLMXBC.EQ.4) GO TO 500
      RETURN
C     INTERNAL WALL AT K=KW (L=LW)
  300 IF(J.GT.JKLW) RETURN
      IF(LORKW.GT.0) L2=LORKW
      GO TO 80
C     SYMMETRY B.C. FOR VELOCITY COMPONENT NORMAL TO UPPER BOUNDARY
  500 DO 600 L=1,LORKMX
  600 F(L,KLMXBC)=0.0
      RETURN
      END
*DECK DZY
      SUBROUTINE DZY(KORL)
*CALL RHSBCC
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL COUNT
      COMMON TURMU(500,60)
      DIMENSION SS(13,40),T(13),U(40),V(40),W(40),E(40)
      EQUIVALENCE (SS,A),(U,B),(V,B(1,2,1)),(W,B(1,3,1)),(E,B(1,4,1))
     1,(T(1),T1,B(1,5,1)),(T(2),T2),(T(3),T3),(T(4),T4)
     2,(T(5),T5),(T(6),T6),(T(7),T7),(T(8),T8),(T(9),T9),(T(10),T10)
     3,(T(11),T11),(T(12),T12),(T(13),T13)
C     VISCOUS STRESSES( INNER FIRST DERIVATIVES) FOR VISCOUS CROSS
C     TERMS IN A K,L PLANE
      IF(IW.EQ.2) GO TO 300
      K=KORL
      CALL ZZM(J,K,1,LMAX)
      DO 100 L=1,LMAX
      CALL YYM(J,L,K,K)
      DO 100 N=1,4
  100 XX(L,N)=YY(K,N)
      DO 200 L=1,LMAX
      DO 200 N=1,4
  200 YY(L,N)=XX(L,N)
      GO TO 600
  300 L=KORL
      CALL YYM(J,L,1,KMAX)
      DO 400 K=1,KMAX
      CALL ZZM(J,K,L,L)
      DO 400 N=1,4
  400 XX(K,N)=ZZ(L,N)
      DO 500 K=1,KMAX
      DO 500 N=1,4
      ZZ(K,N)=YY(K,N)
  500 YY(K,N)=XX(K,N)
  600 K=KORL
      QW=0.0
      DO 800 L=1,LORKMX
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      RA = 1./Q(KL,1,J)
      Q2 = (Q(KL,2,J)**2 + Q(KL,3,J)**2 + Q(KL,4,J)**2)*RA
      TT= GD*(Q(KL,5,J) - 0.5*Q2)*RA
      CALL VISCOF(TT,RMUE)
      VNU = RMUE+TURMU(KL,J)
      GKAP = RMUE+PRTR*TURMU(KL,J)
      RJ = 1./Q(KL,6,J)
      G1=RJ*VNU
      DO 700 N=1,3
      G2=YY(L,N)*G1
      T(N)=ZZ(L,1)*G2
      T(N+3)=ZZ(L,2)*G2
  700 T(N+6)=ZZ(L,3)*G2
      U(L) = Q(KL,2,J)*RA
      V(L) = Q(KL,3,J)*RA
      W(L) = Q(KL,4,J)*RA
      E(L) = Q(KL,5,J)*RA-.5*(U(L)**2+V(L)**2+W(L)**2)
      SS(1,L)=T1+T5+T9
      SS(2,L)=SS(1,L)+T1/3.0
      SS(3,L)=T2-2./3.*T4
      SS(4,L)=T3-2./3.*T7
      SS(5,L)=SS(1,L)+T5/3.0
      SS(6,L)=T4-2./3.*T2
      SS(7,L)=T6-2./3.*T8
      SS(8,L)=SS(1,L)+T9/3.
      SS(9,L)=T7-2./3.*T3
      SS(10,L)=T8-2./3.*T6
      SS(11,L)=SS(6,L)*V(L)+SS(9,L)*W(L)
      SS(12,L)=SS(3,L)*U(L)+SS(10,L)*W(L)
      SS(13,L)=SS(4,L)*U(L)+SS(7,L)*V(L)
      SS(1,L)=GKPR*GKAP*SS(1,L)/VNU
  800 CONTINUE
      I=I+1
      IF(I.EQ.4) I=1
C     LATERAL LIMIT OF INTERNAL AND LOWER WALLS
      KLWP=KORLW
      IF(KORLW.EQ.0) KLWP=KORLMX
      DO 1000 L=1,LORKMX
C     CENTRAL DIFFERENCES AT INTERIOR POINTS
      NDIFAC=1
      L1=L+1
      L2=L-1
C     ONE-SIDED DIFFERENCES AT BOUNDARY POINTS
C     LOWER BOUNDARY L (K) =1
      IF(L.NE.1) GO TO 801
      L2=L
      NDIFAC=2
      GO TO 890
C     UPPER BOUNDARY LMAX (KMAX)
  801 IF(L.NE.LORKMX) GO TO 805
      L1=L
      NDIFAC=2
      GO TO 890
C     INTERNAL WALLS
C     TEST FOR INTERNAL WALL NORMAL TO L (K) DIRECTION THAT WOULD DISRUPT
C     CENTRAL DIFFERENCES IN THAT DIRECTION
  805 IF(LORKW.EQ.0) GO TO 890
      LKWP=LORKW+1
      IF(L.NE.LORKW .AND.  L.NE.LKWP) GO TO 890
C     TEST FOR INTERNAL WALL PARALLEL TO L (K) DIRECTION
      IF(KORLW) 850,850,810
C     ONE-SIDED DIFFERENCES AT UPPER EDGE OF THE INTERNAL WALL
C     PARALLEL TO L (K) DIRECTION
  810 IF(L.NE.LKWP) GO TO 850
      IF(J.GT.JKLW .OR. K.NE.KORLW .OR. K.NE.KORLW+1) GO TO 850
      L2=L
      NDIFAC=2
      GO TO 890
C     ONE-SIDED DIFFERENCES ON FACES OF THE INTERNAL WALL NORMAL TO THE
C     L (K) DIRECTION
  850 IF(K.GT.KLWP.OR. J.GT.JLKW) GO TO 890
      IF(L.EQ.LORKW) L1=L
      IF(L.EQ.LKWP) L2=L
      NDIFAC=2
  890 CONTINUE
      DO 900 N=1,13
  900 T(N)=SS(N,L)
C     DERIVATIVES
      DU = U(L1)-U(L2)
      DV = V(L1)-V(L2)
      DW = W(L1)-W(L2)
      DEI = E(L1)-E(L2)
      DU2=(U(L1)**2-U(L2)**2)*0.5
      DV2=(V(L1)**2-V(L2)**2)*0.5
      DW2=(W(L1)**2-W(L2)**2)*0.5
C     BOUNDARY POINT DERIVATIVES
      IF(L.NE.1) GO TO 930
C     LOWER BOUNDARY L (K)=1
C     TEST B.C. OPTIONS
      IF(LK1BC) 980,980,910
C      TEST FOR LOWER WALL
  910 IF(K.LE.KLWP .AND. J.GE.JLK1WL .AND. J.LE.JLK1WU) GO TO 990
C     SYMMETRY CONDITIONS OUTSIDE WALL REGION
      IF(IW-1) 915,915,920
C     SYMMETRY PLANE Z=CONSTANT
  915 DWFAC=1.0
      DVFAC=0.0
      GO TO 970
C     SYMMETRY PLANE Y=CONSTANT
  920 DWFAC=0.0
      DVFAC=1.
      GO TO 970
C     UPPER BOUNDARY LMAX (KMAX)
  930 IF(L.NE.LORKMX) GO TO 990
C     TEST B.C. OPTIONS
      LKTEMP=  LKMXBC+ 1
      GO TO (980,990,980,920,915,980),LKTEMP
C     INTERNAL WALLS ARE TAKEN CARE OF BY ONE-SIDED DIFFERENCES AS
C     SET UP PREVIOUSLY
C     SYMMETRY CONDITIONS ON DERIVATIVES
  970 DU=0.
      DU2=0.
      DV2=0.
      DW2=0.
      DEI=0.
      DV=DV*DVFAC
      DW=DW*DWFAC
      GO TO 990
C     SET VISCOUS STRESSES TO ZERO
  980 DO 985 N=1,5
      C(L,N,I)=0.0
  985 CONTINUE
      GO TO 1000
C     CORRECTION FACTOR FOR ONE-SIDED DIFFERENCES
  990 IF(NDIFAC.NE.2) GO TO 995
      DU=DU*2.0
      DV=DV*2.0
      DW=DW*2.0
      DU2=DU2*2.0
      DV2=DV2*2.0
      DW2=DW2*2.0
      DEI=DEI*2.0
  995 CONTINUE
C     COMPUTE VISCOUS STRESSES
      C(L,1,I)=0.0
      C(L,2,I)=T2*DU+T3*DV+T4*DW
      C(L,3,I)=T6*DU+T5*DV+T7*DW
      C(L,4,I)=T9*DU+T10*DV+T8*DW
      C(L,5,I)=T1*DEI+T2*DU2+T5*DV2+T8*DW2+T11*DU+T12*DV+T13*DW
 1000 CONTINUE
      RETURN
      END
*DECK DZYDYZ
      SUBROUTINE DZYDYZ(IWW)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL RHSBCC
C     VISCOUS RIGHT HAND SIDE TERMS THAT CONTAIN CROSS DERIVATIVES
C     IN THE K,L OR L,K DIRECTIONS
      DRE = HD/(RE*DZ1*DY1*2.0)
      IW=IWW
      CALL SW
      I=0
      KLXP=KORLMX+1
      DO 400 J=JB,JMAX
      KLWT=KORLW
      IF(KORLW.EQ.0) KLWT=-5
      DO 300 K=1,KLXP
      IF(K.LT.KLXP) CALL DZY(K)
      IF(K.EQ.1) GO TO 300
      IF(K.EQ.2) CALL DYZBC1(1)
      IF(K.GT.2.AND.K.LT.KLXP) CALL DYZ
      IF(K.EQ.KLWT+2) CALL DYZBC1(2)
      IF(K.EQ.KLXP) CALL DYZBCX(1)
      IF(K.EQ.KLWT+1) CALL DYZBCX(2)
      DO 200 L=1,LORKMX
      IF(IW.EQ.1) KL=(L-1)*ND+K-1
      IF(IW.EQ.2) KL=(K-2)*ND+L
      DO 200 N = 2,5
  200 S(KL,N,J) = S(KL,N,J)+F(L,N)*DRE
  300 CONTINUE
  400 CONTINUE
      RETURN
      END
*DECK DZZDYY
      SUBROUTINE DZZDYY(IWW)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL COUNT
*CALL RHSBCC
      COMMON TURMU(500,60)
      COMMON/VISC/S0(40),S1(40),S2(40),S3(40),S4(40),S5(40),S6(40),U(40)
     1  ,V(40),W(40),E(40),RR(40)
C     VISCOUS RIGHT HAND SIDE TERMS THAT CONTAIN ONLY UNIDIRECTIONAL
C     DERIVATIVES IN EITHER THE K OR L DIRECTION.
      IW=IWW
      CALL SW
      LKWR=LORKW+1
      KLWR=KORLW+1
      GKPR = GAMMA/PR
      PRTR = PR/0.9
      DRE = HD/(RE*DZORDY**2)
      DO 2000 J=JB,JMAX
      DO 2000 K = KLORLL,KUORLU
      INTWAL=0
      IF(LORKW.GT.0.AND.J.LE.JLKW) INTWAL=1
      IF(KORLW.GT.0.AND.K.GE.KLWR) INTWAL=0
      IF(IW.EQ.2) GO TO 700
      CALL ZZM(J,K,1,LMAX)
      GO TO 900
  700 CALL YYM(J,K,1,KMAX)
      DO 800 N=1,4
      DO 800 L=1,KMAX
  800 ZZ(L,N)=YY(L,N)
C     SPECIFY BOUNDARY CONDITION ON DIMENSIONLESS WALL HEAT TRANS-
C   FER RATE QW, NORMALIZED BY RHOINF*CINF**3.
  900 QW=0.0
      QW=QW*SQRT(ZZ(1,1)**2 + ZZ(1,2)**2 + ZZ(1,3)**2)
      DO 1000 L = 1,LORKMX
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      RA = 1./Q(KL,1,J)
      Q2 = (Q(KL,2,J)**2 + Q(KL,3,J)**2 + Q(KL,4,J)**2)*RA
      TT = GD*(Q(KL,5,J) - 0.5*Q2)*RA
      CALL VISCOF(TT,RMUE)
      VNU = RMUE+TURMU(KL,J)
      GKAP = RMUE+PRTR*TURMU(KL,J)
      RJ = 1./Q(KL,6,J)
      S0(L) = (ZZ(L,1)**2+ZZ(L,2)**2+ZZ(L,3)**2)*RJ
      S1(L) = (S0(L)+ZZ(L,1)**2/3.*RJ)*VNU
      S2(L) = (S0(L)+ZZ(L,2)**2/3.*RJ)*VNU
      S3(L) = (S0(L)+ZZ(L,3)**2/3.*RJ)*VNU
      S4(L) = (ZZ(L,1)*ZZ(L,2)/3.*RJ)*VNU
      S5(L) = (ZZ(L,1)*ZZ(L,3)/3.*RJ)*VNU
      S6(L) = (ZZ(L,2)*ZZ(L,3)/3.*RJ)*VNU
      S0(L) = S0(L)*GKPR*GKAP
      U(L) = Q(KL,2,J)*RA
      V(L) = Q(KL,3,J)*RA
      W(L) = Q(KL,4,J)*RA
      E(L) = Q(KL,5,J)*RA-.5*(U(L)**2+V(L)**2+W(L)**2)
 1000 CONTINUE
      R2=0.
      R3=0.
      R4=0.
      R5=0.
      DO 1100 L = 1,LMORKM
      L1 = L+1
      T1 = S1(L1)+S1(L)
      T2 = S2(L1)+S2(L)
      T3 = S3(L1)+S3(L)
      T4 = S4(L1)+S4(L)
      T5 = S5(L1)+S5(L)
      T6 = S6(L1)+S6(L)
      T16 = S0(L1)+S0(L)
      DU = U(L1)-U(L)
      DV = V(L1)-V(L)
      DW = W(L1)-W(L)
      DEI = E(L1)-E(L)
      F2 = T1*DU+T4*DV+T5*DW
      F3 = T4*DU+T2*DV+T6*DW
      F4 = T5*DU+T6*DV+T3*DW
C    CONSERVATIVE DIFFERENCING OF VISCOUS STRESS WORK TERMS IN ENERGY EQ
      DU2=U(L1)**2-U(L)**2
      DV2=V(L1)**2-V(L)**2
      DW2=W(L1)**2-W(L)**2
      DUV=U(L1)*V(L1)-U(L)*V(L)
      DUW=U(L1)*W(L1)-U(L)*W(L)
      DVW=V(L1)*W(L1)-V(L)*W(L)
      F5=T16*DEI+(T1*DU2+T2*DV2+T3*DW2)*0.5+T4*DUV+T5*DUW+T6*DVW
      F(L,1) = 0.
      F(L,2) = F2-R2
      F(L,3) = F3-R3
      F(L,4) = F4-R4
      F(L,5) = F5-R5
C     VISCOUS TERMS IN ENERGY EQUATION AT LOWER SURFACE OF INTERMEDIATE WALL
      IF(INTWAL.EQ.0) GO TO 1040
      IF(L.EQ.LORKW) F(L,5)=-2.0*(R5+DZORDY*RE*QW)
 1040 CONTINUE
      R2 = F2
      R3 = F3
      R4 = F4
      R5 = F5
C     VISCOUS TERMS IN ENERGY EQ UPPER SURFACES INTERMEDIATE WALL
      IF(INTWAL.EQ.0) GO TO 1050
      IF(L.EQ.LKWR) F(L,5)=2.0*(R5+DZORDY*RE*QW)
 1050 CONTINUE
      IF(L.GT.1) GO TO 1100
      R21=R2
      R31=R3
      R41=R4
      R51=R5
 1100 CONTINUE
C     DEFINE VISCOUS RHS TERMS AT THE LOWER BOUNDARY L=1 FOR
C     LK1BC=1 (EITHER WALL OR SYMMETRY BC)
      IF(LK1BC.EQ.0) GO TO 1300
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1 SURFACE
      IF(KORLW) 1110,1110,1120
C     TEST FOR WALL AT L=1
 1110 IF(JLK1WL.LE.J.AND.J.LE.JLK1WU) GO TO 1150
C     SYMMETRY
      GO TO 1140
C     TEST FOR WALL AT L=1
 1120 IF(K.LE.KORLW.AND.(JLK1WL.LE.J.AND.J.LE.JLK1WU)) GO TO 1150
C     VISCOUS TERMS AT SYMMETRY PLANE PORTION OF LOWER BDY L=1
 1140 F(1,1)=0.0
      F(1,2)=2.0*R21
      F(1,3)=2.0*R31
      F(1,4)=2.0*R41
      F(1,5)=2.0*R51
      F(1,5-IW)=0.0
      GO TO 1300
C     ADIABATIC WALL BC AT WALL PORTION OF LOWER BDY L=1
 1150 DO 1200 N=1,4
 1200 F(1,N)=0.0
      F(1,5)=2.0*(R51+DZORDY*RE*QW)
 1300 CONTINUE
C     DEFINE VISCOUS RHS TERMS AT THE UPPER BOUNDARY L=LORKMX
C     FOR LKMXBC=1(WALL),2(FREESTREAM), OR 5(OUTFLOW).
      DO 1400 N=1,5
 1400 F(LORKMX,N)=0.0
      IF(LKMXBC.EQ.5) GO TO 1500
      IF(LKMXBC-3) 1450,1650,1650
 1450 IF(LKMXBC-1) 1700,1600,1500
 1500 F(LORKMX,2)=-2.0*R2
      F(LORKMX,3)=-2.0*R3
      F(LORKMX,4)=-2.0*R4
      F(LORKMX,5)=-2.0*R5
      GO TO 1700
 1600 F(LORKMX,5)=-2.0*(R5+DZORDY*RE*QW)
      GO TO 1700
C     DEFINE VISCOUS RHS TERMS AT THE UPPER BOUNDARY L=LORKMX WHEN THE
C     LATTER IS A SYMMETRY PLANE(LKMXBC=3 OR 4).
 1650 F(LORKMX,1)=0.0
      F(LORKMX,2)=-2.0*F2
      F(LORKMX,3)=-2.0*F3
      F(LORKMX,4)=-2.0*F4
      F(LORKMX,5)=-2.0*F5
      F(LORKMX,LKMXBC)=0.0
 1700 CONTINUE
      DO 1800 L = LLORKL,LUORKU
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      DO 1800 N = 1,5
 1800 S(KL,N,J) = S(KL,N,J)+F(L,N)*DRE
 2000 CONTINUE
      RETURN
      END
*DECK EIGEN
      SUBROUTINE EIGEN
*CALL BASE
*CALL VARS
C
C    COMPUTE EIGENVALUES
C
C  IF CNBR INPUTED GT 0.0 & THEN DT IS SET FOR THAT COURANT NUMBER
C  IF CNBR INPUTED EQ 0.0 &THE INPUTED DT IS USED AND THE MAXIMUM
C  COURANT NUMBER IS FOUND
C
      JT = 1
      KT = 1
      LT = 1
      SIGMAX = 0.
      DO 15 L=LL,LU
      DO 15 K=KAL,KU
      KL = (L-1)*ND+K
      CALL XXM(K,L,1,JMAX)
      DO 15 J=JB,JMAX
      CALL YYM(J,L,K,K)
      CALL ZZM(J,K,L,L)
      RR = 1./Q(KL,1,J)
      U = Q(KL,2,J)*RR
      V = Q(KL,3,J)*RR
      W = Q(KL,4,J)*RR
      SNDSP = SQRT(GD*(Q(KL,5,J)*RR-.5*(U*U+V*V+W*W)))
      SIGA = ABS(XX(J,4)+U*XX(J,1)+V*XX(J,2)+W*XX(J,3))+SNDSP*SQRT(
     1  XX(J,1)**2+XX(J,2)**2+XX(J,3)**2)
      SIGB = ABS(YY(K,4)+U*YY(K,1)+V*YY(K,2)+W*YY(K,3))+SNDSP*SQRT(
     1  YY(K,1)**2+YY(K,2)**2+YY(K,3)**2)
      SIGC = ABS(ZZ(L,4)+U*ZZ(L,1)+V*ZZ(L,2)+W*ZZ(L,3))+SNDSP*SQRT(
     1  ZZ(L,1)**2+ZZ(L,2)**2+ZZ(L,3)**2)
      SIGABC = AMAX1(SIGA,SIGB,SIGC)
      SIGABC=SIGABC/Q(KL,6,J)
      IF (SIGABC-SIGMAX) 15,15,10
   10 JT = J
      KT = K
      LT = L
      SIGMAX = SIGABC
   15 CONTINUE
      DXM = AMIN1(DX1,DY1,DZ1)
      IF (CNBR.GT.0.0) GO TO 20
      GO TO 30
   20 DT = CNBR*DXM/SIGMAX
      DTFAC=1.0
      HD = .5*DT
      HDX = HD/DX1
      HDY = HD/DY1
      HDZ = HD/DZ1
      WRITE (6,25)CNBR,JT,KT,LT,DT
   25 FORMAT(1H0,3HCN=,E20.10,3X,11H,AT J,K,L =,3I5,3X,25HTHE TIME STEP
     1IS RESET TO,E20.10)
   30 CONTINUE
      CNBR = DT*SIGMAX/DXM
      WRITE (6,35)JT,KT,LT,CNBR,SMU
   35 FORMAT(1H ,30HJ,K,L, AND MAX COURANT NBR    ,3I5,E14.5,2X,4HSMU=
     1    1PE12.5)
      RETURN
      END
*DECK FILTRX
      SUBROUTINE FILTRX(K,L,J1,J2)
*CALL BASE
*CALL VARS
*CALL BTRID
      DIMENSION E(5,5)
      KL = (L-1)*ND+K
      JA = J1+1
      JBB = J2-1
      CALL XXM(K,L,J1,J2)
      DO 15 J=J1,J2
      R1 = XX(J,1)*HDX
      R2 = XX(J,2)*HDX
      R3 = XX(J,3)*HDX
      R4 = XX(J,4)*HDX
      CALL AMATRX(E,J,KL,R1,R2,R3,R4)
      DO 10 M=1,5
      DO 10 N=1,5
   10 D(J,N,M) = E(N,M)
   15 CONTINUE
      DO 30 J=JA,JBB
      RJ = Q(KL,6,J)
      RR = RM * .5*(RJ+Q(KL,6,J-1))
      RF = RM * .5*(RJ+Q(KL,6,J+1))
      DO 25 N=1,5
      DO 20 M=1,5
      A(J,N,M) = -D(J-1,N,M)
      B(J,N,M) = 0.0
   20 C(J,N,M) = D(J+1,N,M)
      A(J,N,N) = A(J,N,N)-RR
      B(J,N,N) = RJ+RF+RR
      C(J,N,N) = C(J,N,N)-RF
   25 F(J,N) = S(KL,N,J)
   30 CONTINUE
      RETURN
      END
*DECK FILTRY
      SUBROUTINE FILTRY(J,L,K1,K2)
*CALL BASE
*CALL VARS
*CALL BTRID
      DIMENSION E(5,5)
      KA = K1+1
      KB = K2-1
      CALL YYM(J,L,K1,K2)
      DO 15 K=K1,K2
      KL = ((L-1)*ND)+K
      R1 = YY(K,1)*HDY
      R2 = YY(K,2)*HDY
      R3 = YY(K,3)*HDY
      R4 = YY(K,4)*HDY
      CALL AMATRX(E,J,KL,R1,R2,R3,R4)
      DO 10 M=1,5
      DO 10 N=1,5
   10 D(K,N,M) = E(N,M)
   15 CONTINUE
      DO 30 K=KA,KB
      KL = ((L-1)*ND)+K
      RJ = Q(KL,6,J)
      RR = RM * .5*(RJ+Q(KL-1,6,J))
      RF = RM * .5*(RJ+Q(KL+1,6,J))
      DO 25 N=1,5
      DO 20 M=1,5
      A(K,N,M) = -D(K-1,N,M)
      B(K,N,M) = 0.0
   20 C(K,N,M) = D(K+1,N,M)
      A(K,N,N) = A(K,N,N)-RR
      B(K,N,N) = RJ+RF+RR
      C(K,N,N) = C(K,N,N)-RF
   25 F(K,N) = S(KL,N,J)
   30 CONTINUE
      RETURN
      END
*DECK FILTRZ
      SUBROUTINE FILTRZ(J,K,L1,L2)
*CALL BASE
*CALL VARS
*CALL BTRID
      DIMENSION E(5,5)
      LA = L1+1
      LB = L2-1
      CALL ZZM(J,K,L1,L2)
      DO 15 L=L1,L2
      KL = (L-1)*ND+K
      R1 = ZZ(L,1)*HDZ
      R2 = ZZ(L,2)*HDZ
      R3 = ZZ(L,3)*HDZ
      R4 = ZZ(L,4)*HDZ
      CALL AMATRX(E,J,KL,R1,R2,R3,R4)
      DO 10 M=1,5
      DO 10 N=1,5
   10 D(L,N,M) = E(N,M)
   15 CONTINUE
      DO 30 L=LA,LB
      KL = (L-1)*ND+K
      RJ = Q(KL,6,J)
      RR = RM * .5*(RJ+Q(KL-ND,6,J))
      RF = RM * .5*(RJ+Q(KL+ND,6,J))
      DO 25 N=1,5
      DO 20 M=1,5
      A(L,N,M) = -D(L-1,N,M)
      B(L,N,M) = 0.0
   20 C(L,N,M) = D(L+1,N,M)
      A(L,N,N) = A(L,N,N)-RR
      B(L,N,N) = RJ+RF+RR
      C(L,N,N) = C(L,N,N)-RF
   25 F(L,N) = S(KL,N,J)
   30 CONTINUE
      RETURN
      END
*DECK FLUXVE
      SUBROUTINE FLUXVE(J,KL,R1,R2,R3,R4)
*CALL BASE
*CALL VARS
      RR = 1./Q(KL,1,J)
      U = Q(KL,2,J)*RR
      V = Q(KL,3,J)*RR
      W = Q(KL,4,J)*RR
      QS = R4+R1*U+R2*V+R3*W
      PP = GAMI*(Q(KL,5,J)-.5*Q(KL,1,J)*(U*U+V*V+W*W))
      FV(1) = Q(KL,1,J)*QS
      FV(2) = Q(KL,2,J)*QS+R1*PP
      FV(3) = Q(KL,3,J)*QS+R2*PP
      FV(4) = Q(KL,4,J)*QS+R3*PP
      FV(5) = (Q(KL,5,J)+PP)*QS-R4*PP
      RETURN
      END
*DECK GDINTF
      SUBROUTINE GDINTF(J,K1,K2,L1,L2,NMX,DINTGL)
C     DOUBLE INTEGRAL OF FHAT OVER ETA(K) AND ZETA(L) BETWEEN
C     LIMITS (K1,K2) AND (L1,L2)
*CALL BASE
*CALL VARS
*CALL BTRID
      DIMENSION DINTGL(5)
      DO 20 L=L1,L2
      KL0=(L-1)*ND
      DO 15 K=K1,K2
      KL=KL0+K
C  METRICS XIX,XIY,XIZ
      M=K
      J1=J
      CALL XXM(M,L,J1,J1)
      CALL FLUXVE(J,KL,XX(J,1),XX(J,2),XX(J,3),0.)
      DO 10 N=1,NMX
   10 A(K,N,1)=FV(N)
C     REDEFINE A ARRAY AS FLUX DEFECT RELATIVE TO FREESTREAM FLUX
C     IN CASES WHERE AN EXTERNAL FLOW IS PRESENT.  FREESTREAM CONDITIONS
C     ARE ASSUMED TO EXIST AT (J,K,L)=(1,1,LMAX),I.E. AT J=1,KL=LMAX,
C     EXCEPT FOR PURELY INTERNAL FLOW FOR WHICH THERE IS A WALL AT
C     L=LMAX OR AT K=KMAX OR BOTH.
C     CALL FLUXVE(1,LMAX,XX(J,1),XX(J,2),XX(J,3),0.)
C     DO 11 N=1,NMX
C  11 A(K,N,1)=A(K,N,1)-FV(N)
   15 CONTINUE
      DO 20 N=1,NMX
      CALL QDRTR(B(L,N,1),A(1,N,1),DY1,K1,K2)
   20 CONTINUE
      DO 30 N=1,NMX
      CALL QDRTR(ENTGRL,B(1,N,1),DZ1,L1,L2)
   30 DINTGL(N)=ENTGRL
      RETURN
      END
*DECK GDINTW
      SUBROUTINE GDINTW(IWW,J1,J2,LORK,KORL1,KORL2,GFS)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL RHSBCC
      DIMENSION GFS(5)
C     DOUBLE INTEGRAL OF TOTAL STRESS OVER A WALL BOUNDARY SURFACE
C     LORK= CONSTANT BETWEEN LIMITS (J1,J2) AND (KORL1,KORL2).
      GKPR=GAMMA/PR
      IW=IWW
      L=LORK
      NDD=ND
      IF(IW.EQ.2) NDD=1
      CALL SW
      IF(L.EQ.LORKW.OR.L.EQ.LORKMX) NDD=-NDD
      DO 300 J=J1,J2
      DO 200 K=KORL1,KORL2
      IF(IW-1) 10,10,20
   10 KL=(L-1)*ND+K
      CALL ZZM(J,K,L,L)
      GO TO 30
   20 KL=(K-1)*ND+L
      CALL YYM(J,K,L,L)
      DO 25 N=1,4
   25 ZZ(L,N)=YY(L,N)
   30 KL1=KL
      KL2=KL+NDD
      KL3=KL2+NDD
      RA=1./Q(KL,1,J)
      Q2=(Q(KL,2,J)**2+Q(KL,3,J)**2+Q(KL,4,J)**2)*RA
      T=GD*(Q(KL,5,J)-0.5*Q2)*RA
      CALL VISCOF(T,RMUE)
      VNU=RMUE
      GKAP=RMUE
C     RECONSTRUCT GEOMETRIC JACOBIAN FROM VOLUME-AVERAGED ONE
      RJ=3./(4.*Q(KL1,6,J)-Q(KL2,6,J))
      S0 = (ZZ(L,1)**2+ZZ(L,2)**2+ZZ(L,3)**2)
      C(K,1,1)=SQRT(S0)
      S0=S0*RJ
      S1 = (S0+ZZ(L,1)**2/3.*RJ)*VNU
      S2 = (S0+ZZ(L,2)**2/3.*RJ)*VNU
      S3 = (S0+ZZ(L,3)**2/3.*RJ)*VNU
      S4 = (ZZ(L,1)*ZZ(L,2)/3.*RJ)*VNU
      S5 = (ZZ(L,1)*ZZ(L,3)/3.*RJ)*VNU
      S6 = (ZZ(L,2)*ZZ(L,3)/3.*RJ)*VNU
      R1=1.0/Q(KL1,1,J)
      R2=1.0/Q(KL2,1,J)
      R3=1.0/Q(KL3,1,J)
      DU=(4.0*Q(KL2,2,J)*R2-3.0*Q(KL1,2,J)*R1-Q(KL3,2,J)*R3)/(2.0*DZ1)
      DV=(4.0*Q(KL2,3,J)*R2-3.0*Q(KL1,3,J)*R1-Q(KL3,3,J)*R3)/(2.0*DZ1)
      DW=(4.0*Q(KL2,4,J)*R2-3.0*Q(KL1,4,J)*R1-Q(KL3,4,J)*R3)/(2.0*DZ1)
      E1=(Q(KL,5,J)-0.5*Q2)*RA
      Q2=R2*(Q(KL2,2,J)**2+Q(KL2,3,J)**2+Q(KL2,4,J)**2)
      E2=(Q(KL2,5,J)-0.5*Q2)*R2
      Q2=R3*(Q(KL3,2,J)**2+Q(KL3,3,J)**2+Q(KL3,4,J)**2)
      E3=(Q(KL3,5,J)-0.5*Q2)*R3
      DE=(4.0*E2-3.0*E1-E3)/(2.0*DZ1)
      GFS(2)=S1*DU+S4*DV+S5*DW
      GFS(3)=S2*DV+S4*DU+S6*DW
      GFS(4)=S3*DW+S5*DU+S6*DV
      GFS(5)=S0*GKPR*GKAP*DE
      CALL FLUXVE(J,KL,ZZ(L,1),ZZ(L,2),ZZ(L,3),0.)
      DO 100 N=2,5
C     CORRECT FOR DIRECTIONALITY OF DU,DV,DW,DE
      IF(L.EQ.LORKW.OR.L.EQ.LORKMX) GFS(N)=-GFS(N)
  100 A(K,N,1)=FV(N)-GFS(N)/RE
  200 CONTINUE
      CALL QDRTR(D(J,1,1),C,DY1,KORL1,KORL2)
      DO 290 N=2,5
      CALL QDRTR(B(J,N,1),A(1,N,1),DY1,KORL1,KORL2)
  290 CONTINUE
  300 CONTINUE
      CALL QDRTR(AREF,D,DX1,J1,J2)
      DO 400 N=2,5
      CALL QDRTR(ENTGRL,B(1,N,1),DX1,J1,J2)
  400 GFS(N)=ENTGRL
      RETURN
      END
*DECK GLOBL
      SUBROUTINE GLOBL(GF)
*CALL BASE
*CALL VARS
*CALL RHSBCC
*CALL BTRID
      COMMON/GLOB/CWFAC,JT,AT
      DIMENSION GF(5),DINTGL(5),DINJMX(5)
C   GLOBL PERFORMANCE PARAMETERS FROM INTEGRAL FORM OF STEADY STATE
C   CONSERVATION LAWS.
C   FREESTREAM CONDITIONS ASSUMED AT GRID POINT (J,K,L)=(1,1,LMAX)
      NMX=5
C     INTEGRAL OVER OUTFLOW BOUNDARY
      J=JMAX
   10 L1=1
      L2=LMAX
      K1=1
      K2=KMAX
      CALL GDINTF(J,K1,K2,L1,L2,NMX,DINTGL)
      DO 15 N=1,NMX
   15 GF(N)=DINTGL(N)
C     CORRECT FOR INTERNAL WALL LW.GT.0
      IF (LW.EQ.0) GO TO 30
      L1=LW
      L2=LW+1
      K1=1
      K2=KMAX
      IF (KW.EQ.0) GO TO 20
      K2=KW+1
      IF (J.LE.JLW) GO TO 20
      IF (J.GT.JKW) GO TO 30
      K1=KW
   20 CALL GDINTF(J,K1,K2,L1,L2,NMX,DINTGL)
      DO 25 N=1,NMX
   25 GF(N)=GF(N)-DINTGL(N)
C     CORRECT FOR INTERNAL WALL KW.GT.0
   30 IF (KW.EQ.0) GO TO 42
      IF (J.GT.JKW) GO TO 42
      K1=KW
      K2=KW+1
      L1=1
      L2=LMAX
      IF (LW.EQ.0) GO TO 35
      L2=LW
   35 CALL GDINTF(J,K1,K2,L1,L2,NMX,DINTGL)
      DO 40 N=1,NMX
   40 GF(N)=GF(N)-DINTGL(N)
   42 CONTINUE
      IF (J.NE.JMAX) GO TO 50
      DO 45 N=1,NMX
   45 DINJMX(N)=GF(N)
   50 IF (J.EQ.1) GO TO 55
      J=1
      GO TO 10
C     NET OUTFLOW INTEGRAL MINUS INFLOW INTEGRAL
   55 DO 60 N=1,NMX
      GF(N)=-DINJMX(N)+GF(N)
   60 CONTINUE
C     INTEGRALS OVER FREESTREAM BOUNDARIES LORKMX WHEN LKMXBC=2 OR 5
      DO 105 IW=1,2
      CALL SW
      L=LORKMX
      IF (LKMXBC.NE.2.AND.LKMXBC.NE.5) GO TO 105
   65 DO 95 J=1,JMAX
      DO 90 K=1,KORLMX
      GO TO (70,75),IW
   70 KL=(L-1)*ND+K
      CALL ZZM(J,K,LMAX,LMAX)
      R1=ZZ(L,1)
      R2=ZZ(L,2)
      R3=ZZ(L,3)
      GO TO 80
   75 KL=(K-1)*ND+L
      CALL YYM(J,K,KMAX,KMAX)
      R1=YY(L,1)
      R2=YY(L,2)
      R3=YY(L,3)
   80 CALL FLUXVE(J,KL,R1,R2,R3,0.)
      DO 85 N=1,NMX
   85 A(K,N,1)=FV(N)
C     SUBTRACT FREESTREAM VALUES TO GET FLUX DEFECT RELATIVE TO
C     FREESTREAM FLUX (SEE COMMENT IN SUBROUTINE GDINTF).
C     IF(LKMXBC.EQ.1  .OR. KLMXBC.EQ.1) GO TO 90
C     CALL FLUXVE(1,LMAX,R1,R2,R3,0.)
C     DO 86 N=1,NMX
C  86 A(K,N,1)=A(K,N,1)-FV(N)
   90 CONTINUE
      DO 95 N=1,NMX
      CALL QDRTR(B(J,N,1),A(1,N,1),1.0,1,KORLMX)
   95 CONTINUE
      DO 100 N=1,NMX
      CALL QDRTR(ENTGRL,B(1,N,1),DX1,1,JMAX)
  100 GF(N)=GF(N)-ENTGRL
  105 CONTINUE
C     SET SWITCH FOR WALL SURFACE INTEGRATION
      IGFSUR=1
      IF(IGFSUR.EQ.0) GO TO 300
      DO 150 N=2,5
      GF(N)=0.
      DINTGL(N)=0.
  150 DINJMX(N)=0.
C     COMPUTE GF(2) TO GF(4) DIRECTLY FROM WALL SURFACE INTEGRATIONS
      K1=1
      K2=KMAX
      IF(KW.GT.0) K2=KW
      IF(L1BC.EQ.1.AND.JL1WU.GT.0) CALL GDINTW(1,JL1WL,JL1WU,1,K1,K2,
     1    DINTGL)
      IF(LMAXBC.EQ.1) CALL GDINTW(1,1,JMAX,LMAX,K1,K2,GF)
      DO 200 N=2,5
  200 GF(N)=GF(N)+DINTGL(N)
      IF(LW.EQ.0) GO TO 220
      CALL GDINTW(1,1,JLW,LW,K1,K2,DINTGL)
      IF(K2.EQ.KW) K2=K2+1
      CALL GDINTW(1,1,JLW,LW+1,K1,K2,DINJMX)
      DO 210 N=2,5
  210 GF(N)=GF(N)+DINTGL(N)+DINJMX(N)
  220 L1=1
      L2=LMAX
      IF(LW.GT.0) L2=LW
      DO 230 N=2,5
      DINTGL(N)=0.
  230 DINJMX(N)=0.
      IF(K1BC.EQ.1.AND.JK1WU.GT.0) CALL GDINTW(2,JK1WL,JK1WU,1,L1,L2,
     1    DINTGL)
      IF(KMAXBC.EQ.1) CALL GDINTW(2,1,JMAX,KMAX,L1,L2,DINJMX)
      DO 240 N=2,5
  240 GF(N)=GF(N)+DINTGL(N)+DINJMX(N)
      IF(KW.EQ.0) GO TO 260
      CALL GDINTW(2,1,JKW,KW,L1,L2,DINTGL)
      IF(L2.EQ.LW) L2=L2+1
      CALL GDINTW(2,1,JKW,KW+1,L1,L2,DINJMX)
      DO 250 N=2,5
  250 GF(N)=GF(N)+DINTGL(N)+DINJMX(N)
  260 DO 270 N=2,5
      GF(N)=-GF(N)
  270 CONTINUE
  300 CONTINUE
C     CHANGE SCALING OF GENERALIZED FORCE VECTOR GF(N)
C     TO REFERENCE VALUES OF INFLOW BOUNDARY MASS FLUX
C     RRHO*RMACH, DYNAMIC PRESSURE .5*RRHO*RMACH**2, AND
C     KINETIC ENGERGY FLUX .5*RRHO*RMACH**3.  REFERENCE AREA
C     IS EQUAL TO SQUARE OF REFERENCE LENGTH.
      AREF=1.0
      FACM=1.0/RMACH
      GF(1)=FACM*GF(1)/AREF
      GF(5)=FACM*GF(5)
      FACM=2.0*FACM*FACM/AREF
      DO 110 N=2,5
  110 GF(N)=FACM*GF(N)
C     NOZZLE DISCHARGE COEFFICEINT CW
C     INTEGRATE MASS FLOW AT THROAT J=JT
      NMX=1
      J=JT
      L1=1
      L2=LMAX
      IF(LW.GT.0) L2=LW
      K1=1
      K2=KMAX
      IF(KW.GT.0) K2=KW
      CALL GDINTF(J,K1,K2,L1,L2,NMX,DINTGL)
      CW=DINTGL(1)/(CWFAC*AT)
      WRITE (6,115)GF,CW
      PRINT 115,GF,CW
  115 FORMAT(4X,3HGF=1P5E11.4,2X,3HCW=1PE11.4)
      RETURN
      END
*DECK INFLTR
      SUBROUTINE INFLTR(ISW,LORK)
*CALL BASE
*CALL BTRID
*CALL RHSBCC
C     FILTER INTERMEDIATE DELTAQ AT INFLOW BOUNDARY
      J=1
      I1=2
      IW=4-ISW
      CALL SW
C     FILTER IN LORK DIRECTION
      I2=LMORKM
      IF(LORKW.GT.0) I2=LORKW-1
      I3=I1-1
      I4=I2+1
      DO 30 N=1,5
      IF(N.EQ.3.OR.N.EQ.4) GO TO 30
      DO 10 I=I1,I2
   10 A(I,1,1)=0.25*(F(I+1,N)+F(I-1,N)+2.0*F(I,N))
      A(1,1,1)=F(1,N)
      A(I2+1,1,1)=F(I2+1,N)
      IF(LK1BC.EQ.1.AND.J.GE.JLK1WL.AND.J.LE.JLK1WU) GO TO 15
      A(1,1,1)=0.5*(F(1,N)+F(2,N))
   15 IF(I2.EQ.LMORKM.AND.(LKMXBC.EQ.3.OR.LKMXBC.EQ.4))
     1  A(I2+1,1,1)=0.5*(F(I2,N)+F(I2+1,N))
      DO 20 I=I3,I4
   20 F(I,N)=A(I,1,1)
   30 CONTINUE
      DO 40 I=I3,I4
      IF(IW.EQ.1) KL=(I-1)*ND+LORK
      IF(IW.EQ.2) KL=(LORK-1)*ND+I
      F(I,3)=VOU(KL)*F(I,2)
   40 F(I,4)=WOU(KL)*F(I,2)
      RETURN
      END
*DECK INITIA
      SUBROUTINE INITIA
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL VISCO
*CALL BTRID
C     COMMON FREESP CONTAINS FREESTREAM PRESSURE AND DENSITY
      COMMON/FREESP/FSP,FSRHO
      COMMON/GLOB/CWFAC,JT,AT
      COMMON/AXISYM/LAXIS
      COMMON/FILTR/INFLT
      COMMON TURMU(500,60)
      PI = 4.*ATAN(1.)
      NC1=0
C
C  READ DATA
C
      READ (5,10)NMAX,JMAX,KMAX,LMAX,LAMIN,INVISC,J1BC,JMAXBC,KPLANE,K1B
     1C,JK1WL,JK1WU,KW,JKW,KMAXBC
      READ (5,10)L1BC,JL1WL,JL1WU,LW,JLW,LMAXBC,NRST,IWRIT,NGRI,NP,KVIS
     1,LVIS,KLVIS,INFLT,ISUTH,NROUT
      READ (5,15)DT,FSMACH,RMACH,RE,PR,RTDEGK,FSP,FST
      READ (5,15)GAMMA,RMUEXP,TW,CNBR,DTFAC,RM,SMU,OMEGA
      READ(5,15) DTMAX
   10 FORMAT(16I5)
   15 FORMAT(8F10.0)
      WRITE (6,20)NMAX,JMAX,KMAX,LMAX,LAMIN,INVISC,J1BC,JMAXBC,KPLANE,K1
     1BC,JK1WL,JK1WU,KW,JKW,KMAXBC
      WRITE (6,25)L1BC,JL1WL,JL1WU,LW,JLW,LMAXBC,NRST,IWRIT,NGRI,NP,KVIS
     1,LVIS,KLVIS,INFLT,ISUTH,NROUT
      WRITE (6,30)DT,FSMACH,RMACH,RE,PR,RTDEGK,FSP,FST
      WRITE (6,35)GAMMA,RMUEXP,TW,CNBR,DTFAC,RM,SMU,OMEGA
      WRITE(6,37) DTMAX
      PRINT    20,NMAX,JMAX,KMAX,LMAX,LAMIN,INVISC,J1BC,JMAXBC,KPLANE,K1
     1BC,JK1WL,JK1WU,KW,JKW,KMAXBC
      PRINT    25,L1BC,JL1WL,JL1WU,LW,JLW,LMAXBC,NRST,IWRIT,NGRI,NP,KVIS
     1,LVIS,KLVIS,INFLT,ISUTH,NROUT
      PRINT    30,DT,FSMACH,RMACH,RE,PR,RTDEGK,FSP,FST
      PRINT    35,GAMMA,RMUEXP,TW,CNBR,DTFAC,RM,SMU,OMEGA
      PRINT     37,DTMAX
   20 FORMAT(122H1    NMAX    JMAX    KMAX    LMAX    LAMIN  INVISC   J1
     1BC   JMAXBC  KPLANE   K1BC    JK1WL   JK1WU    KW      JKW   KMAXB
     1C     /16I8)
   25 FORMAT(//132H     L1BC    JL1WL   JL1WU    LW      JLW   LMAXBC
     1 NRST   IWRIT   NGRI     NP      KVIS    LVIS   KLVIS   INFLT   IS
     2UTH    NROUT         /16I8)
   30 FORMAT(//7X,2HDT,10X,6HFSMACH,10X,5HRMACH
     1,12X,2HRE,13X,2HPR,11X,6HRTDEGK,10X,3HFSP,13X,3HFST,/1P8E15.7)
   35 FORMAT(//5X,5HGAMMA,10X,6HRMUEXP,11X,2HTW,12X,4HCNBR,10X,5HDTFAC,
     113X,2HRM,12X,3HSMU,11X,5HOMEGA/1P8E15.7)
   37 FORMAT(//5X,5HDTMAX/1PE15.7///)
C     RMACH IS REFERENCE MACH NUMBER
C     RE IS REYNOLD'S NUMBER AT REFERENCE CONDITIONS
C      RTDEGK IS REFERENCE TEMPERATURE IN DEGREES KELVIN
C        (USED ONLY IN COMPUTING THE SUTHERLAND VISCOSITY)
C     FSMACH IS FREESTREAM MACH NUMBER
C       FSP IS DIMENSIONLESS FREESTREAM PRESSURE
C       (NORMALIZED BY REFERENCE PRESSURE)
C     FST IS DIMENSIONLESS FREESTREAM TEMPERATURE
C        (NORMALIZED BY REFERENCE TEMPERATURE)
C     INVISC=1 FOR INVISCID FLOW, 0 FOR VISCOUS FLOW
C        (0 IS MANDATORY)
C     LAMIN=1 FOR LAMINAR FLOW, 0 FOR TURBULENT FLOW
C
C     SET BOUNDARY CONDITION OPTIONS
C     KPLANE=1 FOR PLANAR SYMMETRY
C     KPLANE=0 FOR 3D FLOW
C     LMAXBC DEFINES BDY CONDITIONS AT THE SURFACE L=LMAX.  THE BC ARE HELD
C     FIXED AT THEIR INITIAL VALUES FOR LMAXBC=0.
C     LMAXBC=1 SELECTS IMPLICIT VISCOUS WALL BC
C     LMAXBC=2 SELECTS IMPLICIT OUTFLOW BC WITH ENFORCED FREESTREAM
C              PRESSURE,FREESTREAM DENSITY.  THIS OPTION IS RECOMMENDED FOR
C               CASES WHERE THE SURFACE L=LMAX IS A 'FREESTREAM'
C               OUTER COMPUTATIONAL BOUNDARY
C     LMAXBC=3 SELECTS REFLECTIVE SYMMETRY BC SUCH THAT LMAX IS A
C              SYMMETRY PLANE Z=CONSTANT AND THE
C              W-VELOCITY COMPONENT IS ZERO IN THAT SYMMETRY PLANE.
C     LMAXBC=4 SELECTS REFLECTIVE SYMMETRY BC SUCH THAT LMAX IS A
C              SYMMETRY PLANE Y=CONSTANT AND THE V
C              VELOCITY COMPONENT IS ZERO IN THAT SYMMETRY PLANE.
C     LMAXBC=5 SELECTS IMPLICIT OUTFLOW BC IN WHICH ALL VISCOUS STRESSES
C              AND HEAT FLUXES VANISH AT L=LMAX. (NOT RECOMMENDED FOR A
C              FREESTREAM BOUNDARY)
C   NOTE LMAXBC = 0 CAN BE USED TO EMPLOY TIME-LAGGED BDRY CONDITIONS
C   AT LMAX, WHICH CAN BE CODED INTO SUBR. BC.  SIMILARLY TIME-LAGGED BC
C     AT THE OTHER BDYS L=1,K=1,K=KMAX CAN BE CONSTRUCTED WHEN L1BC,K1BC
C     OR KMAXBC ARE ZERO.
C
C     L1BC DEFINES THE BC AT THE SURFACE L=1
C     L1BC=0 MEANS THE BC ARE HELD FIXED AT THEIR INITIAL VALUES.
C     L1BC=1 MEANS IMPLICIT VISCOUS WALL BC ARE APPLIED IN THE WALL
C            REGION(JL1WL.LE.J.LE.JL1WU) AND (1.LE.K.LE.KW OR KMAX).
C     THE REMAINDER OF THE SURFACE L=1 OUTSIDE THIS WALL REGION IS
C     ASSUMED TO COINCIDE WITH THE Z=0 PLANE, AND REFLECTIVE SYMMETRY BC
C      ARE APPLIED SUCH THAT THE W-COMPONENT OF VELOCITY VANISHES IN
C     THE SYMMETRY PLANE.  WHEN KW.GT.0(SEE BELOW), THE DESCRBED WALL
C     REGION EXTENDS OUT ONLY TO K=KW RATHER THAN KMAX.
C
C     KMAXBC DEFINES THE BC AT THE SURFACE K=KMAX.  THE OPTIONS ARE
C            IDENTICAL TO THOSE FOR LMAXBC
C
C     K1BC DEFINES THE BC AT THE SURFACE K=1
C     K1BC=0 MEANS THE BC ARE HELD FIXED AT THEIR INITIAL VALUES.
C     K1BC=1 MEANS IMPLICIT VISCOUS WALL BC ARE APPLIED IN THE WALL
C            REGION(JK1WL.LE.J.LE.JK1WU) AND (1.LE.L.LE.LMAX).
C     THE REMAINDER OF THE SURFACE K=1 OUTSIDE THIS WALL REGION IS
C     ASSUMED TO COINCIDE WITH THE Y=0 PLANE, AND REFLECTIVE SUMMETRY BC
C      ARE APPLIED SUCH THAT THE V-COMPONENT OF VELOCITY VANISHES IN
C     THE SYMMETRY PLANE.  WHEN LW.GT.0(SEE BELOW), THE DESCRIBED WALL
C     REGION EXTENDS OUT ONLY TO L=LW RATHER THAN LMAX.
C
C     TW=0.0 SELECTS IMPLICIT ADIABATIC WALL BC FOR THE ENERGY EQUATION.
C        A NON-ZERO INPUT VALUE OF TW IMPOSES A BOUNDARY CONDITION ON
C        DIMENSIONLESS ABSOLUTE WALL TEMPERATURE EQUAL TO TW.
      ITW=0
      IF(TW.GT.0.0) ITW=1
C
C   SET OPTIONS FOR INFLOW (J=1) AND OUTFLOW(J=JMAX) BOUNDARIES
C   JMAXBC=2 FOR SUBSONIC OUTFLOW WITH FREESTREAM PRESSURE IMPOSED
C   JMAXBC.LT.2 FOR COMPUTED OUTFLOW B.C.
C   J1BC=0 FOR FIXED (LAGGED) B.C. AT J=1.  J1BC.NE.0 USES INITIAL
C   VALUES OF TOTAL PRESSURE, TOTAL ENTHALPY, AND VELOCITY VECTOR
C   DIRECTION COSINES, TOGETHER WITH A CLOSING EQUATION OBTAINED EITHER
C   FROM MOC THEORY(J1BC.LT.0) OR FROM THE CONTINUITY, U-MOMENTUM,OR
C   ENERGY EQUATION(J1BC.GT.0). SEE COMMENTS IN SUBR. BCJMAT
C
C   KVIS = 1(0) SPECIFIES THAT THE VISCOUS TERMS FOR THE ETA(K)
C   DIRECTION ARE TO BE INCLUDED IN (OMITTED FROM) THE COMPUTATION.
C   LVIS = 1(0) SIMILARLY CONTROLS THE COMPUTATION OF THE VISCOUS
C   TERMS FOR THE ZETA(L) COORDINATE DIRECTION
C   KLVIS = 1(0) SIMILARLY CONTROLS THE COMPUTATION OF THE VISCOUS
C   CROSS DERIVATIVE TERMS FOR THE ETA AND ZETA COORDINATE DIRECTIONS.
C
C   SET NORMAL MASS INJECTION RATE WMDOT,NORMALIZED BY (FREESTREAM
C   DENSITY)*(FREESTREAM SOUND SPEED)
      WMDOT =0.0
C   SET FINITE VOL/FINITE-DIFF WEIGHT FACTOR, WTFAC=0.75(1.0) FOR
C   FINITE-VOL (FINITE-DIFF) WEIGHTING OF DELTAQ AT BDY POINTS.
C     WTFAC=0.75 IS MANDATORY
C
      WTFAC=0.75
      INXBC=1
      NC=NC1
      JB=2
      IF(J1BC.NE.0) JB=1
      LL=2
      IF(L1BC.GT.0) LL=1
      LU=LMAX-1
      IF(LMAXBC.GT.0) LU=LMAX
      IF(KPLANE.EQ.1) KMAXBC=0
      IF(KPLANE.EQ.1) K1BC=0
      KAL=2
      IF(K1BC.GT.0) KAL=1
      KU=KMAX-1
      IF(KMAXBC.GT.0) KU=KMAX
      IF(KPLANE.EQ.1) KU=1
      IF(KPLANE.EQ.1) KAL=1
C  COMPUTE CONSTANTS
      ALP=0.
      FSRHO=FSP/FST
      RE = RE/RMACH
C     REDEFINE FSMACH AS DIMENSIONLESS FREESTREAM VELOCITY
C     NORMALIZED BY REFERENCE SOUND SPEED
      FSMACH=FSMACH*SQRT(FST)
      TSUTH=110.33/RTDEGK
      GAMI = GAMMA-1
      GD = GAMMA*GAMI
      HD = .5*DT
      JM = JMAX-1
      KM = KMAX-1
      LM = LMAX-1
      JU=JMAX
      IF(JMAXBC.EQ.0) JU=JM
      ND=KMAX
      ND2 = KMAX*LMAX
      DX1 = 1.
      DY1 = 1.
      DZ1 = 1.
      HDX = HD/DX1
      HDY = HD/DY1
      HDZ = HD/DZ1
      CWFAC=(2./(GAMMA+1.))**(.5*(GAMMA+1.)/(GAMMA-1.))
C     INITIALIZE ARRAYS
      DO 55 J=1,JMAX
      DO 55 KL=1,ND2
      TURMU(KL,J) = 0.0
      DO 50 N=1,5
   50    S(KL,N,J) = 0.0
   55    CONTINUE
C     RESTART SOLUTION START
      READ(2) KMAXR,JMAXR,LMAXR,ITMAXR,LMAXBR,L1BCR,FSMACR,GAMMAR,RER,
     1                                 SMUUR,DTTR,ALPR,CNBRR,PRR
   60 CONTINUE
      READ(2) NC1,TAU,DTT,NK
      READ(2) ((X(KL,J),KL = 1,ND2),J = 1,JMAX),
     1        ((Y(KL,J),KL = 1,ND2),J = 1,JMAX),
     2        ((Z(KL,J),KL = 1,ND2),J = 1,JMAX)
      READ(2)  (((Q(KL,N,J),KL = 1,ND2),N = 1,6),J = 1,JMAX)
      IF (NC1.GT.NRST) GO TO 125
      IF (NC1.NE.NRST) GO TO 60
      WRITE (6,65)KMAXR,JMAXR,LMAXR,ITMAXR,LMAXBR,L1BCR,FSMACR,GAMMAR,RE
     1R,SMUUR,DTTR,ALPR,CNBRR,PRR
   65 FORMAT(5X,57HTHIS IS A RESTART FROM THE FOLLOWING INITIAL CONDITIO
     1NS     /,
     15X,88HKMAXR,JMAXR,LMAXR,ITMAXR,LMAXBR,L1BCR,FSMACR,GAMMAR,RER,SMUR
     1,DTR,ALPR,CNBRR,PRR                /6I10/8E15.8)
C      SWITCH TO IDENTIFY AXIAL SYMMETRY
C
      LAXIS=0
      IF(KPLANE.EQ.0 .AND. ABS(Y(1,1)-Y(KMAX,1)).LT.1.E-6) LAXIS=1
C     FIND NOZZLE THROAT INDEX J=JT WHERE X=0.
      JT=1
      DO 70 J=1,JMAX
      IF (X(1,J).NE.0) GO TO 70
      JT=J
      GO TO 75
   70 CONTINUE
   75 CONTINUE
C     COMPUTE NOZZLE THROAT AREA  -AT
      KT=KMAX
      IF(KW.GT.0) KT=KW
      LT=LMAX
      IF(LW.GT.0) LT=LW
      DO 85 L=1,LT
      KLO=(L-1)*ND
      DO 80 K=1,KT
      KL=KLO+K
      M=K
      CALL XXM(M,L,JT,JT)
      DA=SQRT(XX(JT,1)*XX(JT,1)+XX(JT,2)*XX(JT,2)+XX(JT,3)*XX(JT,3))
   80 A(K,1,1)=DA
      CALL QDRTR(B(L,1,1),A(1,1,1),DY1,1,KT)
   85 CONTINUE
      CALL QDRTR(AT,B(1,1,1),DZ1,1,LT)
      IF(AT.LT.1.E-6) AT=1.
      WRITE (6,90)JT,AT
   90 FORMAT(2X,27HTHROAT LOCATION X=0. AT JT= I3/
     1   2X,15HTHROAT AREA AT=  1PE12.5)
      CALL JACOB
C
C      ADJUST OUTFLOW PLANE VARIABLES TO SATISFY IMPOSED PRESSURE
C      BOUNDARY CONDITION P=FSP WHEN JMAXBC=2
C
      IF(JMAXBC .NE. 2) GO TO 94
      IF(NRST.GT.0) GO TO 94
      J=JMAX
      DO 92 L=1,LMAX
      KLO=(L-1)*ND
      DO 92 K=1,KMAX
      KL=KLO+K
      RV=Q(KL,2,J)**2+Q(KL,3,J)**2+Q(KL,4,J)**2
      POLD=GD*(Q(KL,5,J)-RV/(2.*Q(KL,1,J)))
      PNEW=FSP/POLD
      Q(KL,1,J)=PNEW*Q(KL,1,J)
      Q(KL,5,J)=FSP/GD+RV/(2.*Q(KL,1,J))
   92 CONTINUE
   94 CONTINUE
C
      CALL BC
C
C   COMPUTE HTOT AND PTOT (REFERENCED TO HREF AND PREF, RESP.), AND
C   VELOCITY VECTOR DIRECTION COSINES VOU=V/U, WOU=W/U AT INFLOW
C   BOUNDARY J=1.
C
      IF (J1BC.EQ.0.AND.INXBC.EQ.0) GO TO 115
      J=1
      DO 105 L=1,LMAX
      KL0=(L-1)*ND
      DO 105 K=1,KMAX
      KL=KL0+K
      RR=1.0/Q(KL,1,J)
      XKE=0.5*(Q(KL,2,J)**2 + Q(KL,3,J)**2 + Q(KL,4,J)**2)*RR**2
      ER=Q(KL,5,J)*RR
      HTOT(KL)=GAMI*(GAMMA*ER-GAMI*XKE)
      PTOT(KL)=HTOT(KL)*Q(KL,1,J)*(GAMMA*GAMI*(ER-XKE)/HTOT(KL
     1                                                   ))**(-1.0/GAMI)
      IF (Q(KL,2,J)) 100,95,100
   95 VOU(KL)=1.0/SQRT(RE*FSMACH)
      WOU(KL)=1.0/SQRT(RE*FSMACH)
      VOU(KL)=0.
      WOU(KL)=0.
      GO TO 105
  100 VOU(KL)=Q(KL,3,J)/Q(KL,2,J)
      WOU(KL)=Q(KL,4,J)/Q(KL,2,J)
  105 CONTINUE
      WRITE (6,110)((L,VOU((L-1)*ND+1),WOU((L-1)*ND+1),HTOT((L-1)*ND+1),
     1PTOT((L-1)*ND+1)),L=1,LMAX)
  110 FORMAT(5X,3HK=1/6X,1HL,8X,3HVOU,12X,3HWOU,11X,4HHTOT,
     1  10X,4HPTOT/(5X,I2,1P4E15.7))
  115 CONTINUE
      CALL OUT2(1,1)
      CALL OUT2(1,LMAX)
      IF(KMAX.EQ.1) GO TO 120
      IK=(KMAX-1)/2
      CALL OUT2(IK,1)
      CALL OUT2(IK,LMAX)
      CALL OUT2(KMAX,1)
      CALL OUT2(KMAX,LMAX)
  120 CONTINUE
      RETURN
  125 WRITE (6,130)NC1,NRST
  130 FORMAT(30H  ***ERROR - INPUT VALUE NRST= I2,60H DOES NOT MATCH ANY
     1 RESTART CYCLE NUMBER.  LAST CYCLE READ =   I3)
      STOP
      END
*DECK JACOB
      SUBROUTINE JACOB
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL RHSBCC
      COMMON/AXISYM/LAXIS
      AFAC=.75
      BFAC=.25
C
C  JACOBIAN FORMED
C
      DO 10 J = 1,JMAX
      DO 10 K=1,KMAX
      DO 10 L = 1,LMAX
      KL = (L-1)*ND+K
      CALL DJMET(J,KL,XJ,YJ,ZJ)
      CALL DKMET(J,K,L,KL,XK,YK,ZK)
      CALL DLMET(J,K,L,KL,XL,YL,ZL)
      DINV=XJ*(YK*ZL-ZK*YL)+YJ*(ZK*XL-XK*ZL)+ZJ*(XK*YL-YK*XL)
      Q(KL,6,J)=DINV
   10 CONTINUE
C     UPPER BOUNDARY SUURFACE L=LORKMX
      DO 35 IW=1,2
      CALL SW
      IF(LKMXBC.NE.1) GO TO 35
      L=LORKMX
      DO 30 J=1,JMAX
      DO 30 K=1,KORLMX
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      IF(IW.EQ.1) KL1=KL-ND
      IF(IW.EQ.2) KL1=KL-1
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
   30 CONTINUE
   35 CONTINUE
      IF(LAXIS.EQ.0) GO TO 40
C     AXIAL SYMMETRY AT L=1 FOR ALL K,J
      DO 37 J=1,JMAX
      DO 37 K=1,KMAX
      KL=K
      KL1=K+ND
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
   37 CONTINUE
   40 CONTINUE
C     LOWER BOUNDARY SURFACE LORK=1
      DO 55 IW=1,2
      CALL SW
      IF(JLK1WL.EQ.0) GO TO 55
      IF(JLK1WU.EQ.0) GO TO 55
      L=1
      IF(KORLW.EQ.0) KLIM=KORLMX
      IF(KORLW.GT.0) KLIM=KORLW
      DO 50 J=JLK1WL,JLK1WU
      DO 50 K=1,KLIM
      IF(IW.EQ.1) KL=K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      KL1=KL+ND
      IF(IW.EQ.2) KL1=KL+1
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
   50 CONTINUE
   55 CONTINUE
C     INTERIOR AND EXTERIOR SURFACES OF AN INTERNAL WALL AT L=LORKW
      DO 70 IW=1,2
      CALL SW
      IF(JLKW.EQ.0) GO TO 70
      IF(LORKW.LE.0) GO TO 70
      IF(KORLW.EQ.0) KLIM=KORLMX
      IF(KORLW.GT.0) KLIM=KORLW+1
      L=LORKW
      DO 60 J=1,JLKW
      DO 60 K=1,KLIM
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      IF(IW.EQ.1) KL1=KL-ND
      IF(IW.EQ.2) KL1=KL-1
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
      IF(IW.EQ.1) KL=L*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L+1
      KL1=KL+ND
      IF(IW.EQ.2) KL1=KL+1
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
   60 CONTINUE
   70 CONTINUE
C     EXTERIOR EDGE WHERE TWO INTERNAL WALLS INTERSECT WHEN BOTH WALLS
C     ARE PRESENT AND ONE WALL IS LONGER THAN THE OTHER
      DO 90 IW=1,2
      CALL SW
      IF(KORLW.LE.0) GO TO 90
      IF(LORKW.LE.0) GO TO 90
      IF(JKLW.LE.JLKW) GO TO 90
      L=LORKW+1
      JS=JLKW+1
      KS=KORLW+1
      DO 80 J=JS,JKLW
      DO 80 K=KORLW,KS
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      KL1=KL+ND
      IF(IW.EQ.2) KL1=KL+1
      Q(KL,6,J)=AFAC*Q(KL,6,J)+BFAC*Q(KL1,6,J)
   80 CONTINUE
   90 CONTINUE
      RETURN
      END
*DECK MUTMIX
      SUBROUTINE MUTMIX(L1,L2)
*CALL BASE
*CALL VARS
      COMMON TURMU(500,60)
      COMMON/VISC/TAS(40),UU(40),SNOR(40),TMO(40),TMI(40),TURMT(40),
     1   DS(40),U(40),V(40),W(40),E(40),RR(40)
C     TURBULENT KINEMATIC VISCOSITY FOR MIXING LAYER AND FULLY DEVELOPED
C     JET REGION
C     CALCULATE MAXIMUM VORTICITY AND VELOCITY EXTREMA
      TASMAX=1.E-5
      UMAX=UU(L1)
      UMIN=UU(L1)
      LMIN=L1+1
      DO 1 L=LMIN,L2
      IF(TAS(L).GT.TASMAX) TASMAX=TAS(L)
      IF(UU(L).GT.UMAX) UMAX=UU(L)
      IF(UU(L).LT.UMIN) UMIN=UU(L)
    1 CONTINUE
C     CALCULATE TURBULENT VISCOSITY
      DO 10 L=L1,L2
      IF(TASMAX) 4,4,6
    4 TURMT(L)=0
      GO TO 10
C     MIXING LENGTH
   6  ELMIX=0.13*ABS(UMAX-UMIN)/TASMAX
C     KINEMATIC VISCOSITY OVER RE
      TURMT(L)=(ELMIX**2)*TAS(L)
   10 CONTINUE
      RETURN
      END
*DECK MUTSPF
      SUBROUTINE MUTSPF(IW,LORK,J,KORL,SPFAC)
*CALL BASE
*CALL VARS
C     COMPUTE EXPONENT SPFAC IN EXPONENTIAL DAMPING FACTOR FOR
C     BOUNDARY LAYER ON A WALL AT L(K)=LORK
      GO TO (10,20),IW
   10 KL=((LORK-1)*ND)+KORL
      GO TO 30
   20 KL=((KORL-1)*ND)+LORK
   30 RM1=1./Q(KL,1,J)
      XKE=.5*(Q(KL,2,J)**2+Q(KL,3,J)**2+Q(KL,4,J)**2)*RM1
      TT=GD*(Q(KL,5,J)-XKE)*RM1
      CALL VISCOF(TT,RMUE)
      SPFAC=SQRT(RE*Q(KL,1,J)/RMUE)
      RETURN
      END
*DECK MUTUR
      SUBROUTINE MUTUR
*CALL BASE
*CALL VARS
*CALL RHSBCC
      COMMON TURMU(500,60)
      COMMON/VISC/TAS(40),UU(40),SNOR(40),TMO(40),TMI(40),TURMT(40),
     1   DS(40),U(40),V(40),W(40),E(40),RR(40)
      COMMON/AXISYM/LAXIS
      DATA IFIRST/0/
      IFIRST=IFIRST+1

C
C     CONTROL ROUTINE FOR TURBULENT VISCOSITY TURMU(KL,J)
C
      DO 3000 J=JB,JMAX
      DO 2500 IW=1,2
      IF(KPLANE.EQ.1.AND.IW.EQ.2) GO TO 2500
      IF(LAXIS.EQ.1.AND.KPLANE.EQ.0.AND.IW.EQ.2) GO TO 2500
      CALL SW
      DO 2000 K=KLORLL,KUORLU
   10 DO 20 L=LLORKL,LUORKU
      GO TO (1,2),IW
    1 KL=(L-1)*ND+K
C     ETA METRICS
      CALL YYM(J,L,K,K)
      E(L)=1./SQRT(YY(K,1)**2+YY(K,2)**2+YY(K,3)**2)
      GO TO 3
    2 KL=(K-1)*ND+L
C     ZETA METRICS
      CALL ZZM(J,L,K,K)
      E(L)=1./SQRT(ZZ(K,1)**2+ZZ(K,2)**2+ZZ(K,3)**2)
C     VELOCITY COMPONENTS, TOTAL VELOCITY, AND ELEMENTAL ARC
C     LENGTH ALONG LORK LINE
    3 RHOI=1./Q(KL,1,J)
      U(L)=RHOI*Q(KL,2,J)
      V(L)=RHOI*Q(KL,3,J)
      W(L)=RHOI*Q(KL,4,J)
      UU(L)=SQRT(U(L)**2+V(L)**2+W(L)**2)
      GO TO (6,7),IW
    6 CALL DLMET(J,K,L,KL,XI,YI,ZI)
      GO TO 8
    7 CALL DKMET(J,L,K,KL,XI,YI,ZI)
    8 DS(L)=SQRT(XI**2+YI**2+ZI**2)
   20 CONTINUE
C     VORTICITY COMPONENT NORMAL TO SURFACE ETA=CONSTANT OR ZETA=CONSTANT
      DO 999 L=LLORKL,LUORKU
      GO TO (31,32),IW
   31 KL=(L-1)*ND+K
      D2=.5/DZ1
      GO TO 33
   32 KL=(K-1)*ND+L
      D2=.5/DY1
C     XI METRICS
  33  CALL DJMET(J,KL,XJ,YJ,ZJ)
C     DERIVATIVES OF VELOCITY COMPONENTS IN LORK DIRECTION
      LP=L+1
      LR=L-1
C     INTERMEDIATE WALL IN LORK DIRECTION
      IF(LORKW.LE.0.OR.L.LT.LORKW.OR.L.GT.LORKW+1.OR.J.GT.JLKW)GO TO 100
      IF(KORLW.GT.0.AND.K.GT.KORLW) GO TO 100
      IF(L.EQ.LORKW) GO TO 700
      IF(L.EQ.LORKW+1) GO TO 500
  100 IF(L.NE.1) GO TO 600
      IF(LAXIS.NE.1) GO TO 150
      DU=0.
      DV=V(LP)-V(L)
      DW=W(LP)-W(L)
      GO TO 910
  150 CONTINUE
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1 SURFACE
      IF(KORLW) 200,200,400
C     TEST FOR WALL AT L=1
  200 IF(JLK1WL.LE.J.AND.J.LE.JLK1WU) GO TO 500
C     SYMMETRY
  300 DU=0.0
      IF(IW.EQ.2) GO TO 310
      DV=0.0
      DW=W(LP)-W(L)
      GO TO 910
  310 DV=V(LP)-V(L)
      DW=0.0
      GO TO 910
C     TEST FOR WALL AT L=1
  400 IF(K-KORLW) 200,200,300
C     FORWARD DIFFERENCE
  500 LR=L
      FAC=2.0
      GO TO 900
  600 IF(L.NE.LORKMX) GO TO 800
      IF(LKMXBC.LT.3.OR.LKMXBC.GT.4) GO TO 700
C     SYMMETRY
      DU=0.0
      DV=0.0
      DW=0.0
      LR=LR
      IF(LKMXBC.EQ.4) DW=W(L)-W(LR)
      IF(LKMXBC.EQ.3) DV=V(L)-V(LR)
      GO TO 910
C     BACKWARD DIFFERENCE
  700 LP=L
      FAC=2.0
      GO TO 900
C     CENTRAL DIFFERENCE
  800 FAC=1.0
  900 LR=LR
      LP=LP
      DU=(U(LP)-U(LR))*D2*FAC
      DV=(V(LP)-V(LR))*D2*FAC
      DW=(W(LP)-W(LR))*D2*FAC
  910 CONTINUE
C     VORTICITY COMPONENT
      TAS(L)=ABS(XJ*DU+YJ*DV+ZJ*DW)*E(L)
  999 CONTINUE
C     COMPUTE TURBULENT KINEMATIC VISCOSITY OVER RE,TURMT(L), ALONG THE
C     CURRENT LORK LINE AND STORE TEMPORARILY IN S(KL,IW,J)
C
C     CONFIGURATIONS WITHOUT INTERMEDIATE WALLS (PURE EXTERNAL OR PURE
C     INTERNAL FLOW)
      IF(KW.GT.0 .OR. LW.GT.0) GO TO 1200
C     UPPER WALL PRESENT?
      IF(LKMXBC-1) 1000,1020,1000
C     NO UPPER WALL
 1000 L2=LORKMX
      L1=1
C     LOWER WALL PRESENT?
      IF(LK1BC*JLK1WU) 9000,9000,1001
C     LOWER WALL OR ITS WAKE
 1001 IF(J-JLK1WU) 1002,1002,8000
 1002 IF(J.LT.JLK1WL) GO TO 9000
C
C     WALL BOUNDARY LAYER
C
 7000 CALL MUTSPF(IW,L1,J,K,SPFAC)
      CALL MUTWWK(L1,L2,1,SPFAC)
      GO TO 1990
C
C     WAKE OF A WALL
C
 8000 CALL MUTWWK(L1,L2,0,0.0)
      GO TO 1990
C
C     NO WALLS, WAKES, OR MIXING LAYERS
C
 9000 DO 9100 L=L1,L2
 9100 TURMT(L)=0.0
      GO TO 1990
C
C
C
C     UPPER WALL IS PRESENT
 1020 L1=LORKMX
C     LOWER WALL PRESENT?
      IF(LK1BC*JLK1WU) 1025,1025,1030
C     UPPER WALL BOUNDARY LAYER SPANS ALL LORK
 1025 L2=1
      GO TO 7000
C     BOTH UPPER AND LOWER WALLS ARE PRESENT
 1030 IF(J.LT.JLK1WL) GO TO 1025
C *** THE VARIABLE LEDGE DEFINES THE COMPUTATIONAL SUBDOMAIN
C *** 1.LE.LORK.LT.LEDGE THAT IS SPANNED BY THE LOWER WALL B.L.
C *** OR WAKE AND THE SUBDOMAIN LEDGE.LE.LORK.LE.LORKMX
C *** THAT IS SPANNED BY THE UPPER WALL B.L. LEDGE IS CURRENTLY
C *** TAKEN AT THE POINT OF MAXIMUM TOTAL VELOCITY
      UMAX=0.
      DO 1040 L=1,LORKMX
      IF(UU(L).LT.UMAX) GO TO 1040
      UMAX=UU(L)
      LEDGE=L
 1040 CONTINUE
C     UPPER WALL B.L.
      L2=LEDGE
      CALL MUTSPF(IW,L1,J,K,SPFAC)
      CALL MUTWWK(L1,L2,1,SPFAC)
C     LOWER WALL B.L. OR WAKE
      L1=1
      L2=LEDGE-1
      IF(J-JLK1WU) 7000,7000,8000
C *** LOGIC FOR FULL 3D NOZZLE CONFIGURATION
C *** THAT HAS INTERMEDIATE WALLS AND SIDE PLATES LORKW.GT.0
 1200 CONTINUE
C     DEFINE THE TOLERANCE EPSWAK TO DISTINGUISH WAKES
C     FROM MIXING LAYERS.  IF THE MAXIMUM VELOCITY DEFECT VMIN-VMAX
C     IN A PROFILE SATISFIES VMIN.LT.EPSWAK*VMAX, WE HAVE
C     A WAKE.  OTHERWISE, WE HAVE A MIXING LAYER
      EPSWAK=0.9
C      IS A SIDEWALL/SIDEPLATE PRESENT?
      IF(KORLW) 1300,1300,1210
C     SIDEWALL/SIDEPLATE PRESENT
 1210 IF(K-KORLW) 1300,1220,1225
 1225 IF(K-(KORLW+1))1220,1220,1230
C     REGION OUTSIDE SIDEWALL/SIDEPLATE
C     MIXING LAYER
 1230 CALL MUTMIX(1,LORKMX)
      GO TO 1990
C     IN PLANE OF SIDEWALL/SIDEPLATE K=KORLW,KORLW+1
 1220 IF(J-MAX0(JKLW,JLKW)) 1250,1250,1235
C     REGION DOWNSTREAM OF ALL WALLS
 1235 UMIN=UU(1)
      LMIN=1
      DO 1236 L=2,LORKMX
      IF(UU(L).GT.UMIN) GO TO 1236
      UMIN=UU(L)
      LMIN=L
 1236 CONTINUE
C     WAKE OR MIXING LAYER? TEST UMIN AGAINST UINF
      IF(UMIN-EPSWAK*UU(LORKMX)) 1240,1230,1230
C     WAKE
 1240 CALL MUTWWK(LMIN,LORKMX,0,0.0)
      CALL MUTWWK(LMIN,1,0,0.0)
      GO TO 1990
C     REGION OF SIDEWALL/SIDEPLATE
C     BOUNDARY LAYER ON UPPER SURFACE OF WALL
 1250 CALL MUTSPF(IW,LORKW+1,J,K,SPFAC)
      CALL MUTWWK(LORKW+1,LORKMX,1,SPFAC)
C     BOUNDARY LAYER OR VERTICAL WALL SURFACE BELOW L=LORKW
       IF(J-JKLW) 9000,9000,1260
 1260 CALL MUTSPF(IW,LORKW,J,K,SPFAC)
      CALL MUTWWK(LORKW,1,1,SPFAC)
      GO TO 1990
C     IN REGION BETWEEN SYMMETRY PLANE K=1 AND SIDEWALL/SIDEPLATE K=KORLW
C     LOWER WALL PRESENT AT L=1 ?
 1300 IF(JLK1WU) 1310,1310,1500
C     NO LOWER WALL
 1310 IF(J.GT.JLKW) GO TO 1400
C     BOUNDARY LAYER ON LOWER SURFACE LORKW OF INTERMEDIATE WALL
      CALL MUTSPF(IW,LORKW,J,K,SPFAC)
      CALL MUTWWK(LORKW,1,1,SPFAC)
C     BOUNDARY LAYER ON UPPER SURFACE OF INTERMEDIATE WALL
      CALL MUTSPF(IW,LORKW+1,J,K,SPFAC)
      CALL MUTWWK(LORKW+1,LORKMX,1,SPFAC)
      GO TO 1990
C     INTERMEDIATE WALL WAKE OR MIXING LAYER
 1400 UMIN=UU(1)
      LMIN=1
      UMAX=UU(1)
      DO 1401 L=2,LORKMX
      IF(UU(L).GT.UMIN) GO TO 1401
      UMIN=UU(L)
      LMIN=L
 1401 CONTINUE
      UMAX=AMIN1(UU(1),UU(LORKMX))
C     WAKE OR MIXING LAYER?
      IF(UMIN-EPSWAK*UMAX) 1410,1410,1230
C     WAKE  - IF BROAD ENOUGH
 1410 IF((LMIN.LE.3).OR.(LMIN.GE.(LMAX-2))) GO TO 1230
C     INTERMEDIATE WALL WAKE CENTERED AT LMIN
      CALL MUTWWK(LMIN,1,0,0.0)
      CALL MUTWWK(LMIN,LORKMX,0,0.0)
      GO TO 1990
 1500 CONTINUE
C     LOWER WALL PRESENT AT L=1
      IF(J.GT.MAX0(JLKW,JLK1WU)) GO TO 1700
C     REGION CONTAINING WALLS
      IF(JLKW-JLK1WU) 1520,1520,1510
 1510 WRITE(6,1515)
      PRINT 1515
 1515 FORMAT( 2X,45H*****ERROR FLAG IN SUBROUTINE MUTUR ****
     1,/2X,51H**THIS SUBROUTINE IS NOT PRESENTLY CODED TO COMPUTE
     2,/2X,51H**THE TURBULENT VISCOSITY TURMU FOR PLUGGED NOZZLES
     3,/2X,54H**WHERE THE CENTRAL PLUG AT L(K)=1 IS SHORTER THAN THE
     4,/2X,48H**UPPER NOZZLE WALL AT L(K)=LW(KW).  THE CURRENT
     5,/2X,54H**LOGIC REQUIRES JL1WU .GE. JLW WHEN BOTH ARE POSITIVE
     6,/2X,51H**AND NON-ZERO.  ADDITIONAL LOGIC WOULD BE REQUIRED
     7,/2X,43H**AFTER STATEMENT 1510 TO HANDLE THE REGION
     8,/2X,26H**.0.LE.JL1WU.LE.J.LE.JLW.
     9,/2X,3H***    )
      STOP
 1520 IF(J.GT.JLKW) GO TO 1600
C     REGION J UPSTREAM OF INTERMEDIATE WALL TRAILING EDGE
C     BOUNDARY LAYER ON OUTER SURFACE OF INTERMEDIATE WALL
      CALL MUTSPF(IW,LORKW+1,J,K,SPFAC)
      CALL MUTWWK(LORKW+1,LORKMX,1,SPFAC)
C     BOUNDARY LAYERS INSIDE NOZZLE
      L1=LORKW
      L2=1
C     BOUNDARY LAYER ON UPPER WALL AHEAD OF PLUG
      IF(J.LT.JLK1WL) GO TO 7000
C     BOUNDARY LAYERS ON UPPER WALL AND ON LOWER WALL (PLUG)
      UMAX=UU(1)
      LEDGE=1
      DO 1530 L=2,LORKW
      IF(UU(L).LT.UMAX) GO TO 1530
      UMAX=UU(L)
      LEDGE=L
 1530 CONTINUE
C     UPPER WALL B.L.
      L2=LEDGE
 1540 CALL MUTSPF(IW,L1,J,K,SPFAC)
      CALL MUTWWK(L1,L2,1,SPFAC)
      IF(L1.EQ.1) GO TO 1990
C     B.L. ON LOWER WALL (PLUG)
      L1=1
      GO TO 1540
 1600 CONTINUE
C     SUBREGION BEHIND UPPER NOZZLE WALL AND ABOVE LOWER WALL (PLUG)
      UMAX=UU(1)
      LEDGE=1
      DO 1610 L=2,LORKMX
      IF(UU(L).LT.UMAX) GO TO 1610
      UMAX=UU(L)
      LEDGE=L
 1610 CONTINUE
C B.L. OVER LOWER WALL (PLUG)
      L1=1
      CALL MUTSPF(IW,L1,J,K,SPFAC)
      CALL MUTWWK(L1,LEDGE,1,SPFAC)
C     UPPER WALL WAKE OR MIXING LAYER OUTSIDE PLUG B.L.
      UMIN=UU(LEDGE)
      LMIN=LEDGE
      LEDGE=LEDGE+1
      DO 1620 L=LEDGE,LORKMX
      IF(UU(L).GT.UMIN) GO TO 1620
      UMIN=UU(L)
      LMIN=L
 1620 CONTINUE
C     WAKE OR MIXING LAYER?
      IF(UMIN-EPSWAK*AMIN1(UMAX,UU(LORKMX)))1630,1640,1640
C     UPPER WALL WAKE CENTERED AT LMIN
 1630 CALL MUTWWK(LMIN,LEDGE,0,0.0)
      CALL MUTWWK(LMIN,LORKMX,0,0.0)
      GO TO 1990
C     MIXING LAYER
 1640 CALL MUTMIX(LEDGE,LORKMX)
      GO TO 1990
 1700 CONTINUE
C     REGION DOWNSTREAM OF ALL WALLS
      UMAX=UU(1)
      LEDGE=1
      DO 1710 L=2,LORKMX
      IF(UU(L).LT.UMAX) GO TO 1710
      UMAX=UU(L)
      LEDGE=L
 1710 CONTINUE
      IF(LEDGE-3) 1720,1720,1730
 1720 LEDGE=1
      UMAX=UU(1)
      GO TO 1800
 1730 IF(LEDGE.GE.LORKMX-2) GO TO 1230
C     REGION 1.LE.L.LE.LEDGE.  TEST FOR LOWER WALL (PLUG) WAKE.
      IF(UU(1)-EPSWAK*UMAX) 1740,1750,1750
C     WAKE OF LOWER WALL (PLUG)
 1740 CALL MUTWWK(1,LEDGE,0,0.0)
      GO TO 1800
C     MIXING LAYER IN LOWER REGION
 1750 CALL MUTMIX(1,LEDGE)
 1800 CONTINUE
C     UPPER REGION L.GT.LEDGE. LOOK FOR WAKE OF UPPER NOZZLE WALL.
      UMIN=UU(LEDGE)
      LMIN=LEDGE
      LEDGE=LEDGE+1
      DO 1810 L=LEDGE,LORKMX
      IF(UU(L).GT.UMIN) GO TO 1810
      UMIN=UU(L)
      LMIN=L
 1810 CONTINUE
C     WAKE OR MIXING LAYER?
      IF(LMIN.GT.LEDGE+2) GO TO 1830
C     MIXING LAYER
 1820 CALL MUTMIX(LEDGE,LORKMX)
      GO TO 1990
 1830 IF(UMIN.GE.EPSWAK*AMIN1(UMAX,UU(LORKMX))) GO TO 1820
C     UPPER WALL WAKE CENTERED ABOUT LMIN
      CALL MUTWWK(LMIN,LEDGE,0,0.)
      CALL MUTWWK(LMIN,LORKMX,0,0.)
      GO TO 1990
C     LOAD KINEMATIC VISCOSITY OVER RE TEMPORARILY IN S(KL,IW,1)
 1990 DO 1999 L=LLORKL,LUORKU
      GO TO (1991,1992),IW
 1991 KL=(L-1)*ND+K
      GO TO 1993
 1992 KL=(K-1)*ND+L
 1993 S(KL,IW,1)=TURMT(L)
 1999 CONTINUE
 2000 CONTINUE
 2500 CONTINUE
C     COMBINE INDIVIDUAL VISCOSITY COEFFICIENTS FOR THE K AND L DIRECTIONS
      IF(KPLANE) 2300,2300,2100
 2100 DO 2200 L=LL,LU
      KLO=(L-1)*ND
      DO 2200 K=KAL,KU
      KL=KLO+K
      TURMU(KL,J)=RE*Q(KL,1,J)*S(KL,1,1)
 2200 CONTINUE
      GO TO 3000
 2300 IF(LAXIS.EQ.1) GO TO 2100
      DO 2400 L=LL,LU
      KLO=(L-1)*ND
      DO 2400 K=KAL,KU
      KL=KLO+K
      TURMU(KL,J)=RE*Q(KL,1,J)*SQRT(S(KL,1,1)**2+S(KL,2,1)**2)
      S(KL,1,1)=0.
      S(KL,2,1)=0.
 2400 CONTINUE
 3000 CONTINUE
      RETURN
      END
*DECK MUTWWK
      SUBROUTINE MUTWWK(L1,L2,IWALL,SPFAC)
*CALL BASE
*CALL VARS
      COMMON TURMU(500,60)
      COMMON/VISC/TAS(40),UU(40),SNOR(40),TMO(40),TMI(40),TURMT(40),
     1   DS(40),U(40),V(40),W(40),E(40),FF(40)
      DATA FK,AP,FKK,F27,FCWK,FCKLEB,FELMIX/0.4,26.,.0168,1.6,.5,.3,.26/
C     TURBULENT KINEMATIC VISCOSITY FOR WALL BOUNDARY LAYER (IWALL=1)
C     OR WAKE (IWALL=0) WHERE WALL (OR WAKE CENTER) IS AT L=L1 AND B.L. EDGE
C     OR WAKE EDGE IS AT L=L2
C     ARC LENGTH SNOR, MAXIMUM VORTICITY, VELOCITY EXTREMA, AND
C     VORTICITY FUNCTION FF.
      RA=SPFAC*SQRT(TAS(L1))/AP
      IDL=1
      IF(L1.GT.L2) IDL=-1
      L=L1
      I2=IABS(L2-L1)+1
      IM=1
      SNOR(1)=0.
      TASMAX=TAS(L1)
C  MAX VORTICITY AND ITS LOCATION
      ITASMX=1
      DO 1 I=2,I2
      L=L+IDL
      IF(TAS(L).LE.TASMAX) GO TO 1
      TASMAX=TAS(L)
      ITASMX=I
    1 CONTINUE
C     FIRST MIN. BEYOND THE ABSOLUTE VORTICITY MAX. DETERMINES THE RANGE
C    OF FF(I)
C
C    TO AVOID SPURIOUS PEAKS OUTSIDE THE SENSIBLE BOUNDARY LAYER
C     OR WAKE REGION
C
      IFF= ITASMX
      L=L1
      DO 2 I=2,I2
      L=L+IDL
      IF(I.LE.ITASMX) GO TO 2
      IF(TAS(L).LE.TAS(L-IDL))GO  TO 2
      IFF=I
      GO TO 3
    2 CONTINUE
    3 L=L1
      UMAX=UU(L1)
      UMIN=UU(L1)
      FF(1)=0.
      FFMAX=0.
      SNMAX=.5*(DS(L1)+DS(L1+IDL))
      DO 5 I=2,I2
      L=L+IDL
    5 SNOR(I)=SNOR(I-1)+.5*(DS(L)+DS(L-IDL))
      L=L1
      DO 10 I=2,I2
      L=L+IDL
      IF(UU(L).GT.UMAX) UMAX=UU(L)
      IF(UU(L).LT.UMIN) UMIN=UU(L)
      E1=RA*SNOR(I)
      IF(E1.GT.50) E1=50
      FF(I)=SNOR(I)*TAS(L)*(1.-IWALL*EXP(-E1))
      IF(I.GT.IFF)  GO TO 10
      IF(FF(I).LT.FFMAX) GO TO 9
      IM=I
      FFMAX=FF(I)
      SNMAX=SNOR(I)
      GO TO 10
    9 IF(IM.EQ.1) GO TO 10
      IF(IWALL.NE.1) GO TO 10
C    FOR WALL BOUNDARY LAYER STOP WITH FIRST MAXIMUM FOUND
C
      IF(FF(I).LE.FFMAX) GO TO 15
   10 CONTINUE
   15 CONTINUE
C     INTERPOLATE FOR MAXIMUM POINT OF FF
      IF(IM.EQ.1 .OR. IM.EQ.I2) GO TO 20
      IMP=IM+1
      IMM=IM-1
      BM=.5*(FF(IMP)-FF(IMM))
      AM=BM-FF(IM)
      IF(AM.EQ.0.) GO TO 20
      DIM=.5*BM/AM
      IF(ABS(DIM).GE.1) GO TO 20
      IF(AM.GE.0.) GO TO 20
      FFMAX=FF(IM)-.5*(BM**2)/AM
      BM=.5*(SNOR(IMP)-SNOR(IMM))
      AM=BM-SNOR(IM)
      SNMAX=SNOR(IM)-DIM*(BM-DIM*AM)
   20 CONTINUE
C     OUTER VISCOSITY
      UDIFF=ABS(UMAX-UMIN)
      FWAK=SNMAX*FFMAX
      IF(FFMAX.GT.FCWK*UDIFF) FWAK=SNMAX*((FCWK*UDIFF)**2)/FFMAX
      DO 30 I=1,I2
      FIA=FCKLEB*SNOR(I)/SNMAX
      IF(FIA.GT.1.E5) FIA=1.E5
      TMO(I)=FKK*F27*FWAK/(1.+5.5*FIA**6)
   30 CONTINUE
C     INNER VISCOSITY
      L=L1-IDL
      IF(IWALL.EQ.0) GO TO 50
C     WALL BOUNDARY LAYER
      DO 40 I=1,I2
      L=L+IDL
      E1=RA*SNOR(I)
      IF(E1.GT.50) E1=50
      TMI(I)=TAS(L)*(FK*SNOR(I)*(1.-EXP(-E1)))**2
   40 CONTINUE
      GO TO 70
C     WAKE
   50 DO 60 I=1,I2
      L=L+IDL
      IF(TASMAX) 52,52,54
   52 TMI(I)=0.
      GO TO 60
   54 TMI(I)=TAS(L)*(FELMIX*UDIFF/TASMAX)**2
   60 CONTINUE
C     LOAD VISCOSITIES INTO ARRAY, USING INNER VALUE UNTIL MATCH POINT IS
C     REACHED
   70 I=1
      L=L1
   80 TURMT(L)=TMI(I)
      I=I+1
      L=L+IDL
      IF(I.GT.I2) GO TO 90
      IF(TMI(I).LE.TMO(I)) GO TO 80
   81 TURMT(L)=TMO(I)
      I=I+1
      L=L+IDL
      IF(I.LE.I2) GO TO 81
C      THE TURMT ARRAY CONTAINS THE KINEMATIC VISCOSITY OVER RE
   90 CONTINUE
      RETURN
      END
*DECK OUT2
      SUBROUTINE OUT2(K,L)
*CALL BASE
*CALL VARS
      COMMON/FREESP/FSP,FSRHO
C
C  WRITE VARIABLES AT K,L TAKING OUT METRICS
C
      KL = (L-1)*ND+K
      WRITE (6,10)K,L
   10 FORMAT(1H0,16HVARIABLES AT K= ,I5,8H AND L= ,I5)
      WRITE (6,15)
   15 FORMAT(1H0,4X,1HJ,5X,1HX,10X,1HY,10X,1HZ,7X,8HRHO/RREF,4X,6HU/AREF
     1  ,5X,6HV/AREF,5X,6HW/AREF,5X,6HT/TREF,5X,7HGP/PREF,6X,2HCP)
      WRITE (6,20)
   20 FORMAT(1H0)
      CPS = 1./(.5*GAMMA*FSMACH**2)
      DO 30 J=1,JMAX
      RR = 1./Q(KL,1,J)
      Q1 = Q(KL,1,J)
      Q2 = Q(KL,2,J)*RR
      Q3 = Q(KL,3,J)*RR
      Q4 = Q(KL,4,J)*RR
      Q5 = Q(KL,5,J)
      PP =         GD*(Q5-.5*(Q2**2+Q3**2+Q4**2)*Q1)
      CP = ((PP/FSP)-1.)*CPS
      TT=PP/Q1
      WRITE (6,25)J,X(KL,J),Y(KL,J),Z(KL,J),Q1,Q2,Q3,Q4,TT,PP,CP
   25 FORMAT(1H ,I5,10E11.4)
   30 CONTINUE
      RETURN
      END
*DECK OUTALL
      SUBROUTINE OUTALL(NCIN)
*CALL BASE
*CALL VARS
*CALL COUNT
      COMMON/AXISYM/LAXIS
      KK=KU
      IF(LAXIS.EQ.1) KK=1
      IGLOB=0
      NCOUT=NCIN
      IF(NCIN.GE.0) GO TO 101
      IF(NCIN.LT.0) IGLOB=1
      IF(NCIN.LT.0) NCOUT=-NCOUT
      WRITE(4) NCOUT,TAU,DT,NK
      WRITE(4) ((X(KL,J),KL=1,ND2),J=1,JMAX),
     1         ((Y(KL,J),KL=1,ND2),J=1,JMAX),
     2         ((Z(KL,J),KL=1,ND2),J=1,JMAX)
      WRITE(4) (((Q(KL,N,J),KL=1,ND2),N=1,6),J=1,JMAX)
  101 CONTINUE
      IF(IWRIT.EQ.2) GO TO 40
      WRITE (6,15)NCOUT
   15 FORMAT(//21H  CURRENT FLOW AT NC= I5,28H AS WRITTEN TO RESTART FIL
     1E    /)
      DO 35 J=1,JMAX
      DO 35 K=1,KK
      WRITE (6,20)J,K
   20 FORMAT(1H0,2X,2HJ=,I3,2X,2HK=,I3,2X,1HL,6X,1HX,11X,1HY,11X,1HZ
     1  ,6X,6HR/RREF,5X,6HU/AREF,5X,6HV/AREF,5X,6HW/AREF,5X,6HT/TREF,
     1 5X,6HP/PREF,5X,3HENT)
      DO 30 L=1,LMAX
      KL = (L-1)*ND+K
      R = Q(KL,1,J)
      RR = 1./R
      U = Q(KL,2,J)*RR
      V = Q(KL,3,J)*RR
      W = Q(KL,4,J)*RR
      E = Q(KL,5,J)
      S2=0.0
      IF(ABS(U).GT.1.0E-17) S2=S2+U**2
      IF(ABS(V).GT.1.0E-17) S2=S2+V**2
      IF(ABS(W).GT.1.0E-17) S2=S2+W**2
      PP = GD*(E-.5*R*S2)
      TT=PP*RR
      ENT = PP/(ABS(R))**GAMMA
      WRITE (6,25)L,X(KL,J),Y(KL,J),Z(KL,J),R,U,V,W,TT,PP,ENT
   25 FORMAT(1H ,14X,I3,10F11.6)
   30 CONTINUE
   35 CONTINUE
   40 CONTINUE
      RETURN
      END
*DECK QDRTR
      SUBROUTINE QDRTR(ENTGRL,ENTGRD,DLT,IL,IU)
      DIMENSION ENTGRD(60)
C   1-D TRAPEZOIDAL QUADRATURE. COMPUTES DEFINITE INTEGRAL,ENTGRL,
C   OF INTEGRAND,ENTGRD(I),BETWEEN LIMITS IL,IU.
C
      ENTGRL=0.5*(ENTGRD(IL)+ENTGRD(IU))
       IF((IU-IL).LT.2) RETURN
      ILP=IL+1
      IUM=IU-1
      DO 1 I=ILP,IUM
    1 ENTGRL=ENTGRL+ENTGRD(I)
      ENTGRL=ENTGRL*DLT
      RETURN
      END
*DECK RHS
      SUBROUTINE RHS
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
      IF(LAMIN.EQ.0) CALL MUTUR
      R0 = -HDZ
      LWR=LW+1
      KWR=KW+1
      IF(KW.LE.0) KWR= KMAX+1
      DO 55 J=JB,JMAX
      INTWAL=0
      IF(LW.GT.1.AND.J.LE.JLW) INTWAL=1
      DO 55 K=KAL,KU
      IF(K.GE.KWR) INTWAL=0
      CALL ZZM(J,K,1,LMAX)
      DO 15 L=1,LMAX
      KL = (L-1)*ND+K
      CALL FLUXVE(J,KL,ZZ(L,1),ZZ(L,2),ZZ(L,3),ZZ(L,4))
      DO 10 N=1,5
   10 C(L,N,1) = FV(N)
   15 CONTINUE
      ISYMIN=0
      ISYMAX=0
C     TEST FOR SYMMETRY AT LMAX
      IF(LMAXBC.EQ.3.OR.LMAXBC.EQ.4) ISYMAX=LMAXBC
C     TEST FOR WALL AT L=1
      IF((K.LT.KWR).AND.(JL1WL.LE.J.AND.J.LE.JL1WU)) GO TO 35
C     SYMMETRY
   25 ISYMIN=4
   35 IF (INTWAL.EQ.0) GO TO 40
      CALL DIFFER(R0,1,LW,ISYMIN,0)
      CALL DIFFER(R0,LWR,LMAX,0,ISYMAX)
      GO TO 45
   40 CALL DIFFER(R0,1,LMAX,ISYMIN,ISYMAX)
   45 CONTINUE
      DO 55 L=LL,LU
      KL = (L-1)*ND + K
      DDV= 0.0
      DO 50 N=1,5
   50 S(KL,N,J) = F(L,N) - DDV*Q(KL,N,J)*Q(KL,6,J)
   55 CONTINUE
      IF (KPLANE.EQ.1) GO TO 110
      R0 = -HDY
      LWR=LW+1
      KWR=KW+1
      IF(LW.LE.0) LWR=LMAX+1
      DO 105 J=JB,JMAX
      INTWAL=0
      IF(KW.GT.1.AND.J.LE.JKW) INTWAL=1
      DO 105 L=LL,LU
      IF(L.GE.LWR) INTWAL=0
      CALL YYM(J,L,1,KMAX)
      DO 65 K=1,KMAX
         KL = ((L-1)*ND)+K
         CALL FLUXVE(J,KL,YY(K,1),YY(K,2),YY(K,3),YY(K,4))
      DO 60 N=1,5
   60    C(K,N,1) = FV(N)
   65    CONTINUE
      ISYMIN=0
      ISYMAX=0
C     TEST FOR SYMMETRY AT KMAX
      IF(KMAXBC.EQ.3.OR.KMAXBC.EQ.4) ISYMAX=KMAXBC
C     TEST FOR WALL AT K=1
      IF((L.LT.LWR).AND.(JK1WL.LE.J.AND.J.LE.JK1WU)) GO TO 85
C     SYMMETRY
   75 ISYMIN=3
   85 IF (INTWAL.EQ.0) GO TO 90
      CALL DIFFER(R0,1,KW,ISYMIN,0)
      CALL DIFFER(R0,KWR,KMAX,0,ISYMAX)
      GO TO 95
   90 CALL DIFFER(R0,1,KMAX,ISYMIN,ISYMAX)
   95 DO 100 N=1,5
      DO 100 K=KAL,KU
         KL = ((L-1)*ND)+K
  100    S(KL,N,J) = F(K,N)+S(KL,N,J)
  105    CONTINUE
  110 CONTINUE
      R0 = -HDX
      DO 130 L=LL,LU
      DO 130 K=KAL,KU
      CALL XXM(K,L,1,JMAX)
            KL = ((L-1)*ND)+K
      DO 120 J=1,JMAX
            CALL FLUXVE(J,KL,XX(J,1),XX(J,2),XX(J,3),XX(J,4))
      DO 115 N=1,5
  115       C(J,N,1) = FV(N)
  120       CONTINUE
      CALL DIFFER(R0,1,JMAX,0,0)
      DO 125 J=JB,JMAX
      DO 125 N=1,5
  125       S(KL,N,J) = F(J,N)+S(KL,N,J)
  130       CONTINUE
      IF(INVISC.NE.0) RETURN
C     VISCOUS RIGHT HAND SIDE
C     VISCOUS TERMS FOR L DIRECTION

      IF(LVIS.EQ.1) CALL DZZDYY(1)
      IF(KPLANE.EQ.1) RETURN
C     VISCOUS TERMS FOR K DIRECTION
      IF(KVIS.EQ.1) CALL DZZDYY(2)
C     VISCOUS CROSS DERIVATIVES
      IF(KLVIS.EQ.0) RETURN
      CALL DZYDYZ(1)
      CALL DZYDYZ(2)
      RETURN
      END
*DECK SMOOTH
      SUBROUTINE SMOOTH
*CALL BASE
*CALL VARS
      COMMON/AXISYM/LAXIS
      DIMENSION CDM(5)
C     GLOBALLY FULLY CONSERVATIVE SMOOTHING
      SM0=0.5*SMU/(8.*(3.0-FLOAT(KPLANE)))
      JMM=JM-1

C     XI DIRECTION
      DO 30 L=LL,LU
      DO 30 K=KAL,KU
      KL=((L-1)*ND)+K
      DO 10 N=1,5
   10 CDM(N)=0.0
      SM1=0.0
      DO 15 J=2,JM
   15 SM1=AMAX1(SM1,ABS((Q(KL,6,J+1)-Q(KL,6,J-1))/Q(KL,6,J)))
      SM1=SM0/(1.0+0.25*SM1)
      DO 20 J=2,JMM
      CJ=Q(KL,6,J+1)+Q(KL,6,J)
      DO 20 N=1,5
      CDJN=CJ*(Q(KL,N,J+2)-3.0*(Q(KL,N,J+1)
     1                    -Q(KL,N,J))-Q(KL,N,J-1))
      S(KL,N,J)=S(KL,N,J)-SM1*(CDJN-CDM(N))
      CDM(N)=CDJN
   20 CONTINUE
      DO 25 N=1,5
   25 S(KL,N,JM)=S(KL,N,JM)+SM1*CDM(N)
   30 CONTINUE

C     ZETA(L) DIRECTION OR ETA(K) DIRECTION
      DO 160 IW=1,2
      LORKW=LW
      KORLW=KW
      JLKW=JLW
      KLORLL=KAL
      KUORLU=KU
      LMORKM=LM
      KORLMX=KMAX
      LKMXBC=LMAXBC
      LK1BC=L1BC
      JLK1WL=JL1WL
      JLK1WU=JL1WU
      NDD=ND
      IF (IW.EQ.1) GO TO 40
      IF (KPLANE.EQ.1) GO TO 160
      JLKW=JKW
      LORKW=KW
      KORLW=LW
      KLORLL=LL
      KUORLU=LU
      LMORKM=KM
      KORLMX=LMAX
      LKMXBC=KMAXBC
      LK1BC=K1BC
      JLK1WL=JK1WL
      JLK1WU=JK1WU
      NDD=1
   40 CONTINUE
      DO 150 J=JB,JMAX
      KLXP=KORLMX
      IF(KORLW.GT.0) KLXP=KORLW
      DO 150 K=KLORLL,KUORLU
      IF(LORKW) 41,41,42
   41 M2=1
      L1=2
      L2=LMORKM
      GO TO 45
   42 IF(K.GT.KLXP.OR.J.GE.JLKW) GO TO 41
      M2=2
      L1=2
      L2=LORKW-1
   45 DO 150 M=1,M2
      IF (M-2) 55,50,55
   50 L1=LORKW+2
      L2=LMORKM
      NSYML=0
      GO TO 60
   55 NSYML=0
C  TEST FOR LOWER SYMMETRY PLANE
      IF(LK1BC.EQ.1.AND.(J.LT.JLK1WL.OR.J.GT.JLK1WU.OR.K.GT.KLXP))
     1   NSYML=5-IW
C  TEST FOR UPPER SYMMETRY PLANE
   60 NSYMU=0
      IF(L2.EQ.LMORKM.AND.(LKMXBC.EQ.3.OR.LKMXBC.EQ.4)) NSYMU = LKMXBC
      DO 65 N=1,5
   65 CDM(N)=0.0
      SM1=0.0
      DO 70 L=L1,L2
      KL=((L-1)*ND)+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      K1=KL+NDD
      K3=KL-NDD
   70 SM1=AMAX1(SM1,ABS((Q(K1,6,J)-Q(K3,6,J))/Q(KL,6,J)))
      SM1=SM0/(1.0+0.25*SM1)
      L2M=L2-1
      IF (LAXIS.NE.1) GO TO 75
      IF (L1.EQ.2.AND.IW.EQ.1) GO TO 80
   75 IF (L1.GT.2.OR.NSYML.EQ.0) GO TO 115
C  SMOOTHING AT LOWER SYMMETRY PLANE OR SYMMETRY AXIS
   80 KL=K
      IF(IW.EQ.2) KL=(K-1)*ND+1
      K1=KL+NDD
      K2=K1+NDD
      CJ=Q(K1,6,J) + Q(KL,6,J)
      DO 110 N=1,5
      IF (LAXIS) 95,95,85
   85 IF (IW.NE.1) GO TO 95
      IF (L1-2) 95,90,95
   90 IF (N.EQ.3.OR.N.EQ.4) GO TO 100
      GO TO 105
   95 IF (N-NSYML) 105,100,105
  100 S(KL,N,J)=S(KL,N,J)-6.0*SM1*CJ*Q(KL,N,J)
      CDM(N)=CJ*(Q(K2,N,J)+Q(KL,N,J)+2.*(Q(KL,N,J)
     1  -Q(K1,N,J)))
      GO TO 110
  105 CDM(N)=CJ*(Q(K2,N,J)-Q(K1,N,J)+3.*(Q(KL,N,J)
     1  -Q(K1,N,J)))
      S(KL,N,J)=S(KL,N,J)-2.0*SM1*CDM(N)
  110 CONTINUE
  115 DO 120 L=L1,L2M
      KL=((L-1)*ND)+K
      IF(IW.EQ.2) KL=((K-1)*ND)+L
      K1=KL+NDD
      K2=KL+2*NDD
      K3=KL-NDD
      CJ=Q(K1,6,J)+Q(KL,6,J)
      DO 120 N=1,5
      CDKLN=CJ*(Q(K2,N,J)-3.0*(Q(K1,N,J)
     1                        -Q(KL,N,J))-Q(K3,N,J))
      S(KL,N,J)=S(KL,N,J)-SM1*(CDKLN-CDM(N))
      CDM(N)=CDKLN
  120 CONTINUE
      KL=L2M*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L2
      K1=KL+NDD
      K3=KL-NDD
      CJ=Q(K1,6,J)+Q(KL,6,J)
      DO 145 N=1,5
      IF (NSYMU) 130,125,130
  125 S(KL,N,J)=S(KL,N,J)+SM1*CDM(N)
      GO TO 145
C  SMOOTHING AT UPPER SYMMETRY PLANE
  130 IF (N-NSYMU) 140,135,140
  135 S(K1,N,J)=S(K1,N,J)-6.0*SM1*CJ*Q(K1,N,J)
      CDKLN=CJ*(-(Q(K1,N,J)+Q(K3,N,J))+2.*(Q(KL,N,J)
     1  -Q(K1,N,J)))
      S(KL,N,J)=S(KL,N,J)-SM1*(CDKLN-CDM(N))
      GO TO 145
  140 CDKLN=CJ*(Q(KL,N,J)-Q(K3,N,J)+3.*(Q(KL,N,J)
     1  -Q(K1,N,J)))
      S(KL,N,J)=S(KL,N,J)-SM1*(CDKLN-CDM(N))
      S(K1,N,J)=S(K1,N,J)+2.0*SM1*CDKLN
  145 CONTINUE
  150 CONTINUE
C     IF(IW.NE.1) GO TO 1000
  160 CONTINUE
      RETURN
      END
*DECK STEP
      SUBROUTINE STEP
*CALL BASE
*CALL VARS
*CALL COUNT
*CALL BTRID
      COMMON/FILTR/INFLT
C
   10 FORMAT(1H0,5X,4H NC=,I3,6H TIME=,F10.4,4H DT=,E11.4,
     1 1X,17HMAXIMUM DELTAQ/Q=,1PE10.3,11HAT J,K,L,N=I2,3(1H,,I2))
C
      CALL RHS
      CALL SMOOTH
C
C  COMPUTE L2 RESIDUAL
C
      VOL=0.
      RESM=0.
      RESID = 0.0
      K1=2
      K2=KM
      JMH=JM-1
      KMH=KM-1
      LMH=LM-1
      IF (KPLANE.EQ.0) GO TO 15
      K1=KAL
      K2=KU
      KMH=1
   15 DO 25 J=2,JM
      DO 25 N=1,5
      DO 25 L=2,LM
      DO 25 K=K1,K2
      KL = ((L-1)*ND)+K
      DVOL=1.0/Q(KL,6,J)
      VOL=VOL+DVOL
      RES=ABS(S(KL,N,J))
      IF (RES-RESM) 25,25,20
   20 RESM=RES
      JRESM=J
      KRESM=K
      LRESM=L
      NRESM=N
   25 RESID = RESID+RES**2
      RESID = SQRT(RESID)/(DT+.00005)/(JMH*KMH*LMH*5.)
      NCM=NC-1
      RESOUT = RESID
      ITEROUT = NC
      WRITE (6,30)NCM,RESID,RESM,JRESM,KRESM,LRESM,NRESM
      PRINT 30,NCM,RESID,RESM,JRESM,KRESM,LRESM,NRESM
   30 FORMAT(1H0,2HN=,I4,1X,9HL2 RESID=,1PE10.3,11H MAX RESID=,
     1  1PE10.3,9H J,K,L,N=,I2,3(1H,,I2))
      IF (RESID.LT.1.) GO TO 40
      WRITE (6,35)
   35 FORMAT(41H ****STOPPING**RESIDUAL IS GREATER THAN 1  )
      CALL OUTALL(NC)
      GO TO 120
   40 CONTINUE
C     J SWEEP ALONG (K,L) LINES
      DO 50 L=LL,LU
      DO 50 K=KAL,KU
      KL = ((L-1)*ND)+K
      CALL FILTRX(K,L,1,JMAX)
      IF(J1BC.NE.0) CALL BCXIN(KL)
      CALL BCXOUT(KL)
      IF(J1BC.NE.0) CALL BCJMAT(1,K,L,1)
      IF(JMAXBC.EQ.2) CALL BCJMAT(JMAX,K,L,1)
      CALL BCALJ(K,L)
C
C  S MUST BE ZERO ON B.C.
C
      CALL BTRI(JB,JMAX)
      DO 45 J=JB,JMAX
      DO 45 N=1,5
   45 S(KL,N,J) = F(J,N)*Q(KL,6,J)
   50 CONTINUE
C
C     K SWEEP ALONG (J,L) LINES
C
      IF (KPLANE.EQ.1) GO TO 80
      DO 75 J=JB,JMAX
      DO 75 L=LL,LU
         CALL FILTRY(J,L,1,KMAX)
      CALL BCLK1(2,J,L)
      CALL BCLKMX(2,J,L)
      IF(INVISC.EQ.0) CALL VISMAT(2,J,L)
      IF (J.EQ.1.AND.J1BC.NE.0) GO TO 55
      IF (J.EQ.JMAX.AND.JMAXBC.EQ.2) GO TO 55
      GO TO 65
   55 DO 60 K=KAL,KU
   60 CALL BCJMAT(J,K,L,2)
   65 CONTINUE
      CALL BCALKL(J,L,1)
         CALL BTRI(KAL,KU)
      IF(J*J1BC.EQ.1.AND.INFLT.EQ.1) CALL INFLTR(2,L)
      DO 70 N=1,5
      DO 70 K=KAL,KU
         KL = ((L-1)*ND)+K
   70    S(KL,N,J) = F(K,N)*Q(KL,6,J)
   75    CONTINUE
   80 CONTINUE
C
C
C     L SWEEP ALONG J,K LINES
C
      R0 = -HDZ
      DQ0QMX=0.
      RUMAX=0.
      DO 115 J=JB,JU
C     SET UP TO AVOID DQOQ FOR N=2 IN SPEARATED REGIONS
      IF(KPLANE.EQ.1) GO TO 84
      DO 83 L=LL,LU
      DO 83 K = KAL,KU
      KL=((L-1)*ND)+K
      IF(Q(KL,2,J).LE.RUMAX) GO TO 83
      RUMAX=Q(KL,2,J)
  83  CONTINUE
   84 CONTINUE
      DO 115 K=KAL,KU
            CALL FILTRZ(J,K,1,LMAX)
      CALL BCLK1(1,J,K)
      CALL BCLKMX(1,J,K)
      IF(INVISC.EQ.0) CALL VISMAT(1,J,K)
      IF (J.EQ.1.AND.J1BC.NE.0) GO TO 85
      IF (J.EQ.JMAX.AND.JMAXBC.EQ.2) GO TO 85
      GO TO 95
   85 DO 90 L=LL,LU
   90 CALL BCJMAT(J,K,L,3)
   95 CONTINUE
      CALL BCALKL(J,K,2)
         CALL BTRI(LL,LU)
      IF(J*J1BC.EQ.1.AND.INFLT.EQ.1) CALL INFLTR(3,K)
      DO 110 L=LL,LU
            KL = (L-1)*ND+K
      DO 105 N=1,5
      IF (N.EQ.3.OR.N.EQ.4) GO TO 100
      IF(N.EQ.2.AND.(ABS(Q(KL,2,J)).LT.0.01*RUMAX))GO TO 100
      IF (Q(KL,N,J).EQ.0.) GO TO 100
        DQ0Q=F(L,N)/Q(KL,N,J)
      IF (ABS(DQ0Q).LE.ABS(DQ0QMX)) GO TO 100
        DQ0QMX=DQ0Q
        JDQ=J
        KDQ=K
        LDQ=L
        NDQ=N
  100 CONTINUE
  105 Q(KL,N,J) = F(L,N) + Q(KL,N,J)
  110 CONTINUE
  115       CONTINUE
C
      CALL BC
C
      TAU = TAU + DT
      WRITE (6,10)NC,TAU,DT,DQ0QMX,JDQ,KDQ,LDQ,NDQ
      PRINT 10,NC,TAU,DT,DQ0QMX,JDQ,KDQ,LDQ,NDQ
      RETURN
C     ERROR STOP WHEN RESIDUAL IS GREATER THAN 1
C
  120 STOP
      END
*DECK SW
      SUBROUTINE SW
*CALL RHSBCC
*CALL BASE
      IF(IW.EQ.2) GO TO 100
      LORKMX=LMAX
      KORLMX=KMAX
      LMORKM=LM
      LLORKL=LL
      LUORKU=LU
      KLORLL=KAL
      KUORLU=KU
      LKMXBC=LMAXBC
      KLMXBC=KMAXBC
      LK1BC=L1BC
      KL1BC=K1BC
      LORKW=LW
      KORLW=KW
      JLKW=JLW
      JKLW=JKW
      JLK1WL=JL1WL
      JLK1WU=JL1WU
      JKL1WL=JK1WL
      JKL1WU=JK1WU
      DZORDY=DZ1
      RETURN
  100 LORKMX=KMAX
      KORLMX=LMAX
      LMORKM=KM
      LLORKL=KAL
      LUORKU=KU
      KLORLL=LL
      KUORLU=LU
      LKMXBC=KMAXBC
      KLMXBC=LMAXBC
      LK1BC=K1BC
      KL1BC=L1BC
      KORLW=LW
      LORKW=KW
      JLKW=JKW
      JKLW=JLW
      JKL1WL=JL1WL
      JKL1WU=JL1WU
      JLK1WL=JK1WL
      JLK1WU=JK1WU
      DZORDY=DY1
      RETURN
      END
*DECK VISCOF
      SUBROUTINE VISCOF(TEMP,RMUE)
*CALL VISCO
C
C     DIMENSIONLESS LAMINAR VISCOSITY COEFFICIENT RMUE IS EVALUATED
C     FROM SUTHERLAND LAW WHEN ISUTH=1, AND FROM POWER LAW WHEN ISUTH=0
C
      IF(ISUTH-1) 20,10,30
C     SUTHERLAND LAW
   10 RMUE= (1.0+TSUTH)*TEMP*SQRT(TEMP)/(TEMP+TSUTH)
      RETURN
C     POWER LAW
   20 RMUE=TEMP**RMUEXP
      RETURN
   30 CONTINUE
      RETURN
      END
*DECK VISMAT
      SUBROUTINE VISMAT(IWW,JJ,K)
*CALL BASE
*CALL VARS
*CALL BTRID
*CALL COUNT
*CALL RHSBCC
      COMMON/VISC/S0(40),S1(40),S2(40),S3(40),S4(40),S5(40),S6(40),U(40)
     1  ,V(40),W(40),E(40),RR(40)
      COMMON TURMU(500,60)
      DIMENSION DU(40,5,5)
C     JACOBIAN MATRICES FOR VISCOUS TERMS THAT ARE TREATED IMPLICITLY
      IW=IWW
      J=JJ
      CALL SW
      GKPR = GAMMA/PR
      PRTR = PR/0.9
      DRE = -HD/(RE*DZORDY**2)
      IF (IW.EQ.2) GO TO 10
      CALL ZZM(J,K,1,LMAX)
      GO TO 20
   10 CALL YYM(J,K,1,KMAX)
      DO 15 N=1,4
      DO 15 L=1,KMAX
   15 ZZ(L,N)=YY(L,N)
   20 R3 = 1./3.
      DO 25 L=1,LORKMX
      IF(IW.EQ.1) KL=(L-1)*ND+K
      IF(IW.EQ.2) KL=(K-1)*ND+L
      RR(L) = 1./Q(KL,1,J)
      Q2 = (Q(KL,2,J)**2 + Q(KL,3,J)**2 + Q(KL,4,J)**2)*RR(L)
      T  = GD*(Q(KL,5,J) - 0.5*Q2)*RR(L)
      CALL VISCOF(T,RMUE)
      VNU = RMUE+TURMU(KL,J)
      GKAP = RMUE+PRTR*TURMU(KL,J)
      RJ = 1./Q(KL,6,J)
      S0(L) = (ZZ(L,1)**2+ZZ(L,2)**2+ZZ(L,3)**2)*RJ
      S1(L) = (S0(L)+R3*ZZ(L,1)**2*RJ)*VNU
      S2(L) = (S0(L)+R3*ZZ(L,2)**2*RJ)*VNU
      S3(L) = (S0(L)+R3*ZZ(L,3)**2*RJ)*VNU
      S4(L) = R3*ZZ(L,1)*ZZ(L,2)*RJ*VNU
      S5(L) = R3*ZZ(L,1)*ZZ(L,3)*RJ*VNU
      S6(L) = R3*ZZ(L,2)*ZZ(L,3)*RJ*VNU
      S0(L) = S0(L)*GKPR*GKAP
      RR(L) = 1./Q(KL,1,J)
      U(L) = Q(KL,2,J)*RR(L)
      V(L) = Q(KL,3,J)*RR(L)
      W(L) = Q(KL,4,J)*RR(L)
      E(L) = Q(KL,5,J)*RR(L)
   25 CONTINUE
      DO 30 L=2,LORKMX
      LR = L-1
      C0 = S0(L)+S0(LR)
      C1 = S1(L)+S1(LR)
      C2 = S2(L)+S2(LR)
      C3 = S3(L)+S3(LR)
      C4 = S4(L)+S4(LR)
      C5 = S5(L)+S5(LR)
      C6 = S6(L)+S6(LR)
      D(L,2,1) =  -(C1*U(LR)+C4*V(LR)+C5*W(LR))*RR(LR)
      DU(L,2,1)=  -(C1*U(L )+C4*V(L )+C5*W(L ))*RR(L )
      D(L,2,2) =  C1*RR(LR)
      DU(L,2,2) = C1*RR(L)
      D(L,2,3) = C4*RR(LR)
      DU(L,2,3) = C4*RR(L)
      D(L,2,4) = C5*RR(LR)
      DU(L,2,4) = C5*RR(L)
      D(L,2,5) = 0.0
      DU(L,2,5) = 0.0
      D(L,3,1) = -(C2*V(LR)+C4*U(LR)+C6*W(LR))*RR(LR)
      DU(L,3,1)= -(C2*V(L )+C4*U(L )+C6*W(L ))*RR(L )
      D(L,3,2) = C4*RR(LR)
      DU(L,3,2) = C4*RR(L)
      D(L,3,3) = C2*RR(LR)
      DU(L,3,3) = C2*RR(L)
      D(L,3,4) = C6*RR(LR)
      DU(L,3,4) = C6*RR(L)
      D(L,3,5) = 0.0
      DU(L,3,5) = 0.0
      D(L,4,1) = -(C3*W(LR)+C5*U(LR)+C6*V(LR))*RR(LR)
      DU(L,4,1) =-(C3*W(L )+C5*U(L )+C6*V(L ))*RR(L )
      D(L,4,2) = C5*RR(LR)
      DU(L,4,2) = C5*RR(L)
      D(L,4,3) = C6*RR(LR)
      DU(L,4,3) = C6*RR(L)
      D(L,4,4) = C3*RR(LR)
      DU(L,4,4) = C3*RR(L)
      D(L,4,5) = 0.0
      DU(L,4,5) = 0.0
      D(L,5,1) = -((C1-C0)*U(LR)**2+(C2-C0)*V(LR)**2+(C3-C0)*W(LR)**2
     1  +2.*C4*U(LR)*V(LR)+2.*C5*U(LR)*W(LR)+2.*C6*V(LR)*W(LR)+C0*E(LR))
     2  *RR(LR)
      DU(L,5,1) = -((C1-C0)*U(L)**2+(C2-C0)*V(L)**2+(C3-C0)*W(L)**2+
     1  2.*C4*U(L)*V(L)+2.*C5*U(L)*W(L)+2.*C6*V(L)*W(L)+C0*E(L))*RR(L)
      D(L,5,2) = ((C1-C0)*U(LR)+C4*V(LR)+C5*W(LR))*RR(LR)
      DU(L,5,2)= ((C1-C0)*U(L )+C4*V(L )+C5*W(L ))*RR(L )
      D(L,5,3) = ((C2-C0)*V(LR)+C4*U(LR)+C6*W(LR))*RR(LR)
      DU(L,5,3)= ((C2-C0)*V(L )+C4*U(L )+C6*W(L ))*RR(L )
      D(L,5,4) = ((C3-C0)*W(LR)+C5*U(LR)+C6*V(LR))*RR(LR)
      DU(L,5,4) = ((C3-C0)*W(L)+C5*U(L)+C6*V(L))*RR(L)
      D(L,5,5) = C0*RR(LR)
      DU(L,5,5) = C0*RR(L)
   30 CONTINUE
C     TEST FOR INTERMEDIATE WALL
      INTWAL=0
      IF(LORKW.GT.0.AND.J.LE.JLKW) INTWAL=1
      IF(KORLW.GT.0.AND.K.GE.KORLW+1) INTWAL=0
      DO 40 L=2,LMORKM
      LR = L+1
      DO 35 N=2,5
      DO 35 M=1,5
      IF (INTWAL.EQ.1.AND.(L.EQ.LORKW.OR.L.EQ.LORKW+1)) GO TO 35
      A(L,N,M) = A(L,N,M)+DRE*D(L,N,M)
      B(L,N,M) = B(L,N,M)-DRE*(D(LR,N,M)+DU(L,N,M))
      C(L,N,M) = C(L,N,M)+DRE*DU(LR,N,M)
   35 CONTINUE
   40 CONTINUE
C     IMPLICIT VISCOUS BLOCKS FOR INTERMEDIATE WALL
      IF (INTWAL.EQ.0) GO TO 50
      L=LORKW
      LR=LORKW+1
      DO 45 N=2,5
      DO 45 M=1,5
      A(L,N,M)=A(L,N,M)+2.0*DRE*D(L,N,M)
      B(L,N,M)=B(L,N,M)-2.0*DRE*DU(L,N,M)
      C(L,N,M)=0.0
      A(LR,N,M)=0.0
      B(LR,N,M)=B(LR,N,M)-2.0*DRE*D(LR+1,N,M)
   45 C(LR,N,M)=C(LR,N,M)+2.0*DRE*DU(LR+1,N,M)
   50 CONTINUE
C
C     FILL IN VISCOUS PORTION OF IMPLICIT COEFF. MATRIX FOR
C   IMPLICIT OUTFLOW B.C. AT L=LORKMX (LKMXBC=2 OR 5), OR FOR IMPLICIT
C   ADIABATIC WALL B.C. (LKMXBC=1), OR FOR SYMMETRY BC(LKMXBC=3,4).
C
      IF (LKMXBC.EQ.0) GO TO 65
      L=LORKMX
      DO 60 N=2,5
      IF (LKMXBC.LT.3.OR.LKMXBC.GT.4) GO TO 55
      IF (N.EQ.LKMXBC) GO TO 65
   55 DO 60 M=1,5
      A(L,N,M)=A(L,N,M) + (2.0*DRE)*D(L,N,M)
      B(L,N,M)=B(L,N,M) - (2.0*DRE)*DU(L,N,M)
   60 CONTINUE
   65 CONTINUE
C
C     FILL IN VISCOUS PORTION OF IMPLICIT COEFF MATRIX FOR IMPLICIT
C     ADIABATIC WALL BC FOR JLK1WL.LE.J.LE.JLK1WU AND SYMMETRY BC
C     FOR J OUTSIDE THAT INTERVAL
      L=1
      IF (LK1BC.NE.1) GO TO 100
      NSYM=10
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1
      IF (KORLW) 70,70,80
C     TEST FOR WALL AT L=1
   70 IF (JLK1WL.LE.J.AND.J.LE.JLK1WU) GO TO 85
C     SYMMETRY
   75 IF(IW.EQ.1) NSYM=4
      IF(IW.EQ.2) NSYM=3
      GO TO 85
C     TEST FOR WALL AT L=1
   80 IF (K-KORLW) 70,70,75
   85 DO 95 N=2,5
      IF (N.EQ.NSYM) GO TO 95
      DO 90 M=1,5
      B(L,N,M)=B(L,N,M)-(2.0*DRE)*D(L+1,N,M)
      C(L,N,M)=C(L,N,M)+(2.0*DRE)*DU(L+1,N,M)
   90 CONTINUE
   95 CONTINUE
  100 RETURN
      END
*DECK XXM
      SUBROUTINE XXM(M,L,J1,J2)
*CALL BASE
*CALL VARS
*CALL COUNT
C
C  XI METRICS FORMED FOR A K,L LINE IN J
C
C
C  SYMMETRY
C
      K = M
      KL = (L-1)*ND+K
      DO 10 J = J1,J2
      CALL DKMET(J,K,L,KL,XK,YK,ZK)
      CALL DLMET(J,K,L,KL,XL,YL,ZL)
      XX(J,1) = YK*ZL-ZK*YL
      XX(J,2) = ZK*XL-XK*ZL
      XX(J,3) = XK*YL-YK*XL
      XX(J,4) = 0.
   10 CONTINUE
      RETURN
      END
*DECK YYM
      SUBROUTINE YYM(J,L,K1,K2)
*CALL BASE
*CALL VARS
*CALL COUNT
C
C  ETA METRICS FORMED FOR A J,L LINE IN K
C
      DO 10 K = K1,K2
      KL = (L-1)*ND+K
      CALL DJMET(J,KL,XJ,YJ,ZJ)
      CALL DLMET(J,K,L,KL,XL,YL,ZL)
      YY(K,1) = ZJ*YL-YJ*ZL
      YY(K,2) = XJ*ZL-XL*ZJ
      YY(K,3) = YJ*XL-XJ*YL
      YY(K,4) = 0.
   10 CONTINUE
      RETURN
      END
*DECK ZERODQ
      SUBROUTINE ZERODQ(I)
*CALL BASE
*CALL VARS
*CALL BTRID
C     DEFINE IMPLICIT OPERATOR MATRIX ELEMENTS TO YIELD DELTAQ=0
C     AT GRID POINTS WHERE LAGGED BOUNDARY CONDITIONS ARE TO BE USED
      DO 30 N=1,5
      DO 20 M=1,5
      A(I,N,M)=0.
      B(I,N,M)=0.
   20 C(I,N,M)=0.
      B(I,N,N)=1.0
   30 F(I,N)=0.
      RETURN
      END
*DECK ZZM
      SUBROUTINE ZZM(J,M,L1,L2)
*CALL BASE
*CALL VARS
*CALL COUNT
C
C  ZETA METRICS FORMED FOR A J,K LINE IN L
C
C  SYMMETRY
C
      K = M
      DO 10 L = L1,L2
      KL = (L-1)*ND+K
      CALL DKMET(J,K,L,KL,XK,YK,ZK)
      CALL DJMET(J,KL,XJ,YJ,ZJ)
      ZZ(L,1) = YJ*ZK-ZJ*YK
      ZZ(L,2) = XK*ZJ-XJ*ZK
      ZZ(L,3) = XJ*YK-YJ*XK
      ZZ(L,4) = 0.
   10 CONTINUE
      RETURN
      END
